Explore and analyze the continuous chlorophyll data to be included in the NDFS synthesis manuscript. We will attempt to fit multiple models to the data set and use a model selection process to determine the best one.
# Load packages
library(tidyverse)
library(scales)
library(knitr)
library(mgcv)
library(car)
library(gratia)
library(here)
library(conflicted)
# Source functions
source(here("manuscript_synthesis/src/global_functions.R"))
# Declare package conflict preferences
conflicts_prefer(dplyr::filter())
Display current versions of R and packages used for this analysis:
devtools::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.2.3 (2023-03-15 ucrt)
## os Windows 10 x64 (build 19044)
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate English_United States.utf8
## ctype English_United States.utf8
## tz America/Los_Angeles
## date 2023-08-31
## pandoc 3.1.1 @ C:/Program Files/RStudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
## bslib 0.4.2 2022-12-16 [1] CRAN (R 4.2.2)
## cachem 1.0.8 2023-05-01 [1] CRAN (R 4.2.3)
## callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.2)
## car * 3.1-2 2023-03-30 [1] CRAN (R 4.2.3)
## carData * 3.0-5 2022-01-06 [1] CRAN (R 4.2.1)
## cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.3)
## colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.2.2)
## conflicted * 1.2.0 2023-02-01 [1] CRAN (R 4.2.2)
## crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.1)
## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.2.1)
## digest 0.6.31 2022-12-11 [1] CRAN (R 4.2.2)
## dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.2.3)
## ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.1)
## evaluate 0.21 2023-05-05 [1] CRAN (R 4.2.3)
## fansi 1.0.4 2023-01-22 [1] CRAN (R 4.2.2)
## fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.2.2)
## forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.2.2)
## fs * 1.6.2 2023-04-25 [1] CRAN (R 4.2.3)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.1)
## ggplot2 * 3.4.2 2023-04-03 [1] CRAN (R 4.2.3)
## glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.1)
## gratia * 0.8.1 2023-02-02 [1] CRAN (R 4.2.2)
## gtable 0.3.3 2023-03-21 [1] CRAN (R 4.2.3)
## here * 1.0.1 2020-12-13 [1] CRAN (R 4.2.1)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.2.3)
## htmltools 0.5.5 2023-03-23 [1] CRAN (R 4.2.3)
## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.2.3)
## httpuv 1.6.9 2023-02-14 [1] CRAN (R 4.2.2)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.1)
## jsonlite 1.8.4 2022-12-06 [1] CRAN (R 4.2.2)
## knitr * 1.42 2023-01-25 [1] CRAN (R 4.2.2)
## later 1.3.0 2021-08-18 [1] CRAN (R 4.2.1)
## lattice 0.21-8 2023-04-05 [1] CRAN (R 4.2.3)
## lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.2.1)
## lubridate * 1.9.2 2023-02-10 [1] CRAN (R 4.2.2)
## magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.1)
## Matrix 1.5-4 2023-04-04 [1] CRAN (R 4.2.3)
## memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.1)
## mgcv * 1.8-42 2023-03-02 [1] CRAN (R 4.2.2)
## mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.1)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.1)
## mvnfast 0.2.8 2023-02-23 [1] CRAN (R 4.2.2)
## nlme * 3.1-162 2023-01-31 [2] CRAN (R 4.2.3)
## patchwork 1.1.2 2022-08-19 [1] CRAN (R 4.2.1)
## pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.3)
## pkgbuild 1.4.0 2022-11-27 [1] CRAN (R 4.2.2)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.1)
## pkgload 1.3.2 2022-11-16 [1] CRAN (R 4.2.2)
## prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.1)
## processx 3.8.1 2023-04-18 [1] CRAN (R 4.2.3)
## profvis 0.3.7 2020-11-02 [1] CRAN (R 4.2.1)
## promises 1.2.0.1 2021-02-11 [1] CRAN (R 4.2.1)
## ps 1.7.5 2023-04-18 [1] CRAN (R 4.2.3)
## purrr * 1.0.1 2023-01-10 [1] CRAN (R 4.2.2)
## R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.1)
## Rcpp 1.0.10 2023-01-22 [1] CRAN (R 4.2.2)
## readr * 2.1.4 2023-02-10 [1] CRAN (R 4.2.2)
## remotes 2.4.2 2021-11-30 [1] CRAN (R 4.2.1)
## rlang 1.1.1 2023-04-28 [1] CRAN (R 4.2.3)
## rmarkdown 2.21 2023-03-26 [1] CRAN (R 4.2.3)
## rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.2.1)
## rstudioapi 0.14 2022-08-22 [1] CRAN (R 4.2.1)
## sass 0.4.6 2023-05-03 [1] CRAN (R 4.2.3)
## scales * 1.2.1 2022-08-20 [1] CRAN (R 4.2.1)
## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.1)
## shiny 1.7.4 2022-12-15 [1] CRAN (R 4.2.2)
## stringi 1.7.12 2023-01-11 [1] CRAN (R 4.2.2)
## stringr * 1.5.0 2022-12-02 [1] CRAN (R 4.2.2)
## tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.2.3)
## tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.2.2)
## tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.1)
## tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.2.2)
## timechange 0.2.0 2023-01-11 [1] CRAN (R 4.2.2)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.2.3)
## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.1)
## usethis 2.1.6 2022-05-25 [1] CRAN (R 4.2.1)
## utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.2)
## vctrs 0.6.2 2023-04-19 [1] CRAN (R 4.2.3)
## withr 2.5.0 2022-03-03 [1] CRAN (R 4.2.1)
## xfun 0.39 2023-04-20 [1] CRAN (R 4.2.3)
## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.1)
## yaml 2.3.7 2023-01-23 [1] CRAN (R 4.2.2)
##
## [1] C:/R/win-library/4.2
## [2] C:/Program Files/R/R-4.2.3/library
##
## ──────────────────────────────────────────────────────────────────────────────
# Define file path for processed data
fp_data <- here("manuscript_synthesis/data/processed")
# Import daily average chlorophyll and water temperature data
df_chla_wt <- readRDS(file.path(fp_data, "chla_wt_daily_avg_2013-2019.rds"))
# Import daily average flow data
df_flow <- readRDS(file.path(fp_data, "flow_daily_avg_2013-2019.rds"))
# Create a vector for the factor order of StationCode
sta_order <- c(
"RCS",
"RD22",
"I80",
"LIS",
"STTD",
"LIB",
"RVB"
)
# We will use LIS flow data as a proxy for STTD and RD22 flow data as a proxy
# for I-80
df_flow_c <- df_flow %>%
filter(StationCode %in% c("RD22", "LIS")) %>%
mutate(StationCode = case_match(StationCode, "RD22" ~ "I80", "LIS" ~ "STTD")) %>%
bind_rows(df_flow)
# Prepare chlorophyll and flow data for exploration and analysis
df_chla_wt_c1 <- df_chla_wt %>%
# Join flow data to chlorophyll data
left_join(df_flow_c, by = join_by(StationCode, Date)) %>%
mutate(
# Scale and log transform chlorophyll values
Chla_log = log(Chla * 1000),
# Add Year and DOY variables
Year = year(Date),
DOY = yday(Date)
) %>%
# Add Flow Action Periods
ndfa_action_periods() %>%
mutate(
# Apply factor orders to FlowActionPeriod and StationCode
FlowActionPeriod = factor(FlowActionPeriod, levels = c("Before", "During", "After")),
StationCode = factor(StationCode, levels = sta_order),
# Add a column for Year as a factor for the model
Year_fct = factor(Year)
) %>%
arrange(StationCode, Date)
df_chla_wt_c1 %>%
summarize(
min_date = min(Date),
max_date = max(Date),
num_samples = n(),
.by = c(StationCode, Year, FlowActionPeriod)
) %>%
arrange(StationCode, Year, FlowActionPeriod) %>%
kable()
| StationCode | Year | FlowActionPeriod | min_date | max_date | num_samples |
|---|---|---|---|---|---|
| RCS | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| RCS | 2014 | During | 2014-09-09 | 2014-09-23 | 11 |
| RCS | 2014 | After | 2014-09-24 | 2014-09-24 | 1 |
| RCS | 2015 | Before | 2015-08-17 | 2015-08-20 | 4 |
| RCS | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RCS | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RCS | 2016 | Before | 2016-06-23 | 2016-07-13 | 20 |
| RCS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RCS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RCS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RCS | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RCS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RCS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RCS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RCS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RCS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RCS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RCS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| RD22 | 2013 | Before | 2013-08-15 | 2013-08-21 | 7 |
| RD22 | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| RD22 | 2013 | After | 2013-10-03 | 2013-11-04 | 33 |
| RD22 | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| RD22 | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| RD22 | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| RD22 | 2015 | Before | 2015-07-24 | 2015-08-20 | 28 |
| RD22 | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RD22 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RD22 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| RD22 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RD22 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RD22 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RD22 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RD22 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RD22 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RD22 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RD22 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RD22 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RD22 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RD22 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| I80 | 2013 | Before | 2013-08-15 | 2013-08-21 | 7 |
| I80 | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| I80 | 2013 | After | 2013-10-03 | 2013-11-01 | 30 |
| I80 | 2014 | Before | 2014-07-25 | 2014-09-08 | 43 |
| I80 | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| I80 | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| I80 | 2015 | Before | 2015-07-24 | 2015-08-07 | 15 |
| I80 | 2015 | During | 2015-08-27 | 2015-10-01 | 36 |
| I80 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| I80 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| I80 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| I80 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| I80 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| I80 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| I80 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| I80 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| I80 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| I80 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| I80 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| I80 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| I80 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIS | 2013 | Before | 2013-07-07 | 2013-08-21 | 39 |
| LIS | 2013 | During | 2013-08-22 | 2013-10-02 | 38 |
| LIS | 2013 | After | 2013-10-03 | 2013-11-17 | 34 |
| LIS | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| LIS | 2014 | During | 2014-09-09 | 2014-09-23 | 14 |
| LIS | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| LIS | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIS | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| LIS | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIS | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIS | 2017 | During | 2017-08-29 | 2017-09-09 | 12 |
| LIS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| LIS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| LIS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| LIS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| STTD | 2013 | Before | 2013-08-15 | 2013-08-21 | 7 |
| STTD | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| STTD | 2013 | After | 2013-10-03 | 2013-11-04 | 33 |
| STTD | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| STTD | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| STTD | 2014 | After | 2014-09-24 | 2014-11-08 | 29 |
| STTD | 2015 | Before | 2015-07-27 | 2015-08-20 | 25 |
| STTD | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| STTD | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| STTD | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| STTD | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| STTD | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| STTD | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| STTD | 2017 | During | 2017-08-29 | 2017-09-01 | 4 |
| STTD | 2017 | After | 2017-09-20 | 2017-09-25 | 6 |
| STTD | 2018 | Before | 2018-07-20 | 2018-08-27 | 39 |
| STTD | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| STTD | 2018 | After | 2018-09-27 | 2018-10-15 | 19 |
| STTD | 2019 | Before | 2019-07-26 | 2019-08-25 | 31 |
| STTD | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| STTD | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIB | 2013 | Before | 2013-07-16 | 2013-08-21 | 37 |
| LIB | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| LIB | 2013 | After | 2013-10-03 | 2013-11-16 | 45 |
| LIB | 2014 | Before | 2014-07-25 | 2014-09-08 | 46 |
| LIB | 2014 | During | 2014-09-09 | 2014-09-23 | 15 |
| LIB | 2014 | After | 2014-09-24 | 2014-11-05 | 43 |
| LIB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIB | 2015 | During | 2015-08-21 | 2015-10-01 | 38 |
| LIB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIB | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIB | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| LIB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIB | 2018 | Before | 2018-08-14 | 2018-08-27 | 14 |
| LIB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIB | 2019 | Before | 2019-07-11 | 2019-07-30 | 20 |
| LIB | 2019 | During | 2019-09-05 | 2019-09-21 | 17 |
| LIB | 2019 | After | 2019-09-22 | 2019-11-06 | 39 |
| RVB | 2013 | Before | 2013-07-07 | 2013-08-21 | 46 |
| RVB | 2013 | During | 2013-08-22 | 2013-10-02 | 42 |
| RVB | 2013 | After | 2013-10-03 | 2013-11-17 | 46 |
| RVB | 2014 | Before | 2014-07-25 | 2014-09-07 | 45 |
| RVB | 2014 | During | 2014-09-10 | 2014-09-23 | 14 |
| RVB | 2014 | After | 2014-09-24 | 2014-11-08 | 46 |
| RVB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| RVB | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RVB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| RVB | 2016 | Before | 2016-05-29 | 2016-07-13 | 39 |
| RVB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RVB | 2016 | After | 2016-08-02 | 2016-09-16 | 37 |
| RVB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RVB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RVB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RVB | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RVB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RVB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RVB | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RVB | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RVB | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
Looking at the sample counts and date ranges, we’ll only include years 2015-2019 for the analysis. However, there are still a few Station-Year-FlowPeriod combinations with low sample counts within 2015-2019.
df_chla_wt_c2 <- df_chla_wt_c1 %>%
filter(Year %in% 2015:2019) %>%
mutate(Year_fct = fct_drop(Year_fct))
df_chla_wt_c2 %>%
summarize(
min_date = min(Date),
max_date = max(Date),
num_samples = n(),
.by = c(StationCode, Year, FlowActionPeriod)
) %>%
arrange(StationCode, Year, FlowActionPeriod) %>%
kable()
| StationCode | Year | FlowActionPeriod | min_date | max_date | num_samples |
|---|---|---|---|---|---|
| RCS | 2015 | Before | 2015-08-17 | 2015-08-20 | 4 |
| RCS | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RCS | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RCS | 2016 | Before | 2016-06-23 | 2016-07-13 | 20 |
| RCS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RCS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RCS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RCS | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RCS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RCS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RCS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RCS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RCS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RCS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RCS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| RD22 | 2015 | Before | 2015-07-24 | 2015-08-20 | 28 |
| RD22 | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RD22 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RD22 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| RD22 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RD22 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RD22 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RD22 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RD22 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RD22 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RD22 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RD22 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RD22 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RD22 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RD22 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| I80 | 2015 | Before | 2015-07-24 | 2015-08-07 | 15 |
| I80 | 2015 | During | 2015-08-27 | 2015-10-01 | 36 |
| I80 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| I80 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| I80 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| I80 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| I80 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| I80 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| I80 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| I80 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| I80 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| I80 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| I80 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| I80 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| I80 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIS | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIS | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| LIS | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIS | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIS | 2017 | During | 2017-08-29 | 2017-09-09 | 12 |
| LIS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| LIS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| LIS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| LIS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| STTD | 2015 | Before | 2015-07-27 | 2015-08-20 | 25 |
| STTD | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| STTD | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| STTD | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| STTD | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| STTD | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| STTD | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| STTD | 2017 | During | 2017-08-29 | 2017-09-01 | 4 |
| STTD | 2017 | After | 2017-09-20 | 2017-09-25 | 6 |
| STTD | 2018 | Before | 2018-07-20 | 2018-08-27 | 39 |
| STTD | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| STTD | 2018 | After | 2018-09-27 | 2018-10-15 | 19 |
| STTD | 2019 | Before | 2019-07-26 | 2019-08-25 | 31 |
| STTD | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| STTD | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIB | 2015 | During | 2015-08-21 | 2015-10-01 | 38 |
| LIB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIB | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIB | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| LIB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIB | 2018 | Before | 2018-08-14 | 2018-08-27 | 14 |
| LIB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| LIB | 2019 | Before | 2019-07-11 | 2019-07-30 | 20 |
| LIB | 2019 | During | 2019-09-05 | 2019-09-21 | 17 |
| LIB | 2019 | After | 2019-09-22 | 2019-11-06 | 39 |
| RVB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| RVB | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RVB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| RVB | 2016 | Before | 2016-05-29 | 2016-07-13 | 39 |
| RVB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RVB | 2016 | After | 2016-08-02 | 2016-09-16 | 37 |
| RVB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RVB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RVB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RVB | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RVB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RVB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RVB | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RVB | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RVB | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
Before we run our models, we’ll need to remove the NA
values in the Flow variable. When we remove all records where there are
missing daily average flow values, there are even more
Station-Year-FlowPeriod combinations with low sample counts within
2015-2019.
df_chla_wt_q <- df_chla_wt_c2 %>% drop_na(Flow)
df_chla_wt_q %>%
drop_na(Flow) %>%
summarize(
min_date = min(Date),
max_date = max(Date),
num_samples = n(),
.by = c(StationCode, Year, FlowActionPeriod)
) %>%
arrange(StationCode, Year, FlowActionPeriod) %>%
kable()
| StationCode | Year | FlowActionPeriod | min_date | max_date | num_samples |
|---|---|---|---|---|---|
| RCS | 2015 | Before | 2015-08-17 | 2015-08-20 | 4 |
| RCS | 2015 | During | 2015-08-21 | 2015-09-28 | 39 |
| RCS | 2016 | Before | 2016-06-23 | 2016-07-13 | 8 |
| RCS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RCS | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RCS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RCS | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RCS | 2017 | After | 2017-09-19 | 2017-11-03 | 36 |
| RCS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RCS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RCS | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RCS | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RCS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RCS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| RD22 | 2015 | Before | 2015-07-24 | 2015-08-20 | 28 |
| RD22 | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RD22 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| RD22 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| RD22 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| RD22 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| RD22 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RD22 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RD22 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RD22 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RD22 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RD22 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RD22 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| RD22 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RD22 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| I80 | 2015 | Before | 2015-07-24 | 2015-08-07 | 15 |
| I80 | 2015 | During | 2015-08-27 | 2015-10-01 | 36 |
| I80 | 2015 | After | 2015-10-02 | 2015-11-10 | 40 |
| I80 | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| I80 | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| I80 | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| I80 | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| I80 | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| I80 | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| I80 | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| I80 | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| I80 | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| I80 | 2019 | Before | 2019-07-11 | 2019-08-25 | 46 |
| I80 | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| I80 | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIS | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIS | 2015 | During | 2015-08-21 | 2015-10-01 | 37 |
| LIS | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIS | 2016 | Before | 2016-05-29 | 2016-07-13 | 46 |
| LIS | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIS | 2016 | After | 2016-08-02 | 2016-09-16 | 43 |
| LIS | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIS | 2017 | During | 2017-08-29 | 2017-09-09 | 12 |
| LIS | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| LIS | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| LIS | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| LIS | 2018 | After | 2018-09-27 | 2018-11-11 | 42 |
| LIS | 2019 | Before | 2019-07-19 | 2019-08-25 | 34 |
| LIS | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| LIS | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| STTD | 2015 | Before | 2015-07-27 | 2015-08-20 | 25 |
| STTD | 2015 | During | 2015-08-21 | 2015-10-01 | 37 |
| STTD | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| STTD | 2016 | Before | 2016-06-23 | 2016-07-13 | 21 |
| STTD | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| STTD | 2016 | After | 2016-08-02 | 2016-09-16 | 43 |
| STTD | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| STTD | 2017 | During | 2017-08-29 | 2017-09-01 | 4 |
| STTD | 2017 | After | 2017-09-20 | 2017-09-25 | 6 |
| STTD | 2018 | Before | 2018-07-20 | 2018-08-27 | 39 |
| STTD | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| STTD | 2018 | After | 2018-09-27 | 2018-10-11 | 15 |
| STTD | 2019 | Before | 2019-07-26 | 2019-08-25 | 27 |
| STTD | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| STTD | 2019 | After | 2019-09-22 | 2019-11-06 | 46 |
| LIB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| LIB | 2015 | During | 2015-08-21 | 2015-10-01 | 33 |
| LIB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| LIB | 2016 | Before | 2016-05-29 | 2016-07-13 | 36 |
| LIB | 2016 | During | 2016-07-14 | 2016-08-01 | 19 |
| LIB | 2016 | After | 2016-08-02 | 2016-09-16 | 46 |
| LIB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| LIB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| LIB | 2017 | After | 2017-09-19 | 2017-11-02 | 29 |
| LIB | 2018 | Before | 2018-08-14 | 2018-08-27 | 9 |
| LIB | 2018 | During | 2018-08-28 | 2018-09-08 | 12 |
| LIB | 2018 | After | 2018-09-27 | 2018-10-31 | 15 |
| LIB | 2019 | Before | 2019-07-11 | 2019-07-30 | 18 |
| LIB | 2019 | During | 2019-09-07 | 2019-09-21 | 10 |
| LIB | 2019 | After | 2019-09-22 | 2019-10-25 | 3 |
| RVB | 2015 | Before | 2015-07-06 | 2015-08-20 | 46 |
| RVB | 2015 | During | 2015-08-21 | 2015-10-01 | 42 |
| RVB | 2015 | After | 2015-10-02 | 2015-11-16 | 46 |
| RVB | 2016 | Before | 2016-05-29 | 2016-07-13 | 39 |
| RVB | 2016 | During | 2016-07-14 | 2016-08-01 | 4 |
| RVB | 2016 | After | 2016-08-02 | 2016-09-16 | 37 |
| RVB | 2017 | Before | 2017-07-14 | 2017-08-28 | 46 |
| RVB | 2017 | During | 2017-08-29 | 2017-09-18 | 21 |
| RVB | 2017 | After | 2017-09-19 | 2017-11-03 | 46 |
| RVB | 2018 | Before | 2018-07-13 | 2018-08-27 | 46 |
| RVB | 2018 | During | 2018-08-28 | 2018-09-26 | 30 |
| RVB | 2018 | After | 2018-09-27 | 2018-11-11 | 46 |
| RVB | 2019 | Before | 2019-07-11 | 2019-08-25 | 39 |
| RVB | 2019 | During | 2019-08-26 | 2019-09-21 | 27 |
| RVB | 2019 | After | 2019-09-22 | 2019-11-06 | 36 |
df_chla_wt_q %>%
ggplot(aes(x = FlowActionPeriod, y = Chla, fill = FlowActionPeriod)) +
geom_boxplot() +
facet_grid(rows = vars(Year), cols = vars(StationCode)) +
scale_y_log10(labels = trans_format("log10", math_format(10^.x))) +
annotation_logticks(sides = "l") +
theme_bw() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.position = "top"
)
We’ll try running a GAM including all two-way interactions of the categorical predictors - Year, Flow Action Period, Station - and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_cat <- gam(
Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)
Lets look at the model summary and diagnostics:
summary(m_gam_cat)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.92533 0.08173 109.206 < 2e-16 ***
## Year_fct2016 -0.12304 0.09641 -1.276 0.201987
## Year_fct2017 0.43975 0.08230 5.343 9.72e-08 ***
## Year_fct2018 0.15352 0.08074 1.901 0.057351 .
## Year_fct2019 0.06346 0.08131 0.781 0.435138
## FlowActionPeriodDuring -0.01236 0.06708 -0.184 0.853795
## FlowActionPeriodAfter 0.58644 0.07712 7.605 3.67e-14 ***
## StationCodeRD22 0.70108 0.08498 8.250 2.23e-16 ***
## StationCodeI80 0.84468 0.08759 9.644 < 2e-16 ***
## StationCodeLIS -0.43029 0.08360 -5.147 2.80e-07 ***
## StationCodeSTTD -0.68280 0.08732 -7.820 7.00e-15 ***
## StationCodeLIB -1.25628 0.08535 -14.720 < 2e-16 ***
## StationCodeRVB -1.28890 0.08343 -15.449 < 2e-16 ***
## Year_fct2016:FlowActionPeriodDuring -0.11791 0.07122 -1.656 0.097897 .
## Year_fct2017:FlowActionPeriodDuring -0.07107 0.05402 -1.316 0.188334
## Year_fct2018:FlowActionPeriodDuring -0.03788 0.05025 -0.754 0.451023
## Year_fct2019:FlowActionPeriodDuring -0.48857 0.05151 -9.485 < 2e-16 ***
## Year_fct2016:FlowActionPeriodAfter -0.16200 0.07631 -2.123 0.033843 *
## Year_fct2017:FlowActionPeriodAfter -0.26360 0.04869 -5.414 6.59e-08 ***
## Year_fct2018:FlowActionPeriodAfter -0.15729 0.04867 -3.232 0.001241 **
## Year_fct2019:FlowActionPeriodAfter 0.06520 0.04933 1.322 0.186312
## Year_fct2016:StationCodeRD22 -0.27091 0.09536 -2.841 0.004527 **
## Year_fct2017:StationCodeRD22 -0.60584 0.09026 -6.713 2.23e-11 ***
## Year_fct2018:StationCodeRD22 -0.42354 0.08782 -4.823 1.48e-06 ***
## Year_fct2019:StationCodeRD22 -0.21418 0.08845 -2.421 0.015512 *
## Year_fct2016:StationCodeI80 0.35595 0.09680 3.677 0.000240 ***
## Year_fct2017:StationCodeI80 -0.23252 0.09193 -2.529 0.011478 *
## Year_fct2018:StationCodeI80 -0.49323 0.08953 -5.509 3.88e-08 ***
## Year_fct2019:StationCodeI80 0.45168 0.09016 5.010 5.73e-07 ***
## Year_fct2016:StationCodeLIS 0.75033 0.09371 8.007 1.60e-15 ***
## Year_fct2017:StationCodeLIS 0.23373 0.09054 2.581 0.009882 **
## Year_fct2018:StationCodeLIS -0.08130 0.08738 -0.930 0.352235
## Year_fct2019:StationCodeLIS 0.63379 0.08855 7.157 1.00e-12 ***
## Year_fct2016:StationCodeSTTD 0.68029 0.09599 7.087 1.66e-12 ***
## Year_fct2017:StationCodeSTTD -0.06767 0.10052 -0.673 0.500848
## Year_fct2018:StationCodeSTTD -0.61835 0.09139 -6.766 1.55e-11 ***
## Year_fct2019:StationCodeSTTD -0.82420 0.08995 -9.163 < 2e-16 ***
## Year_fct2016:StationCodeLIB 1.17255 0.09423 12.443 < 2e-16 ***
## Year_fct2017:StationCodeLIB -0.27805 0.09103 -3.054 0.002272 **
## Year_fct2018:StationCodeLIB -1.92894 0.10180 -18.949 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.61755 0.10553 -5.852 5.32e-09 ***
## Year_fct2016:StationCodeRVB 0.35867 0.09694 3.700 0.000219 ***
## Year_fct2017:StationCodeRVB -0.55956 0.08929 -6.267 4.14e-10 ***
## Year_fct2018:StationCodeRVB -0.38660 0.08685 -4.452 8.80e-06 ***
## Year_fct2019:StationCodeRVB -0.65441 0.08842 -7.401 1.70e-13 ***
## FlowActionPeriodDuring:StationCodeRD22 -0.41717 0.06339 -6.581 5.40e-11 ***
## FlowActionPeriodAfter:StationCodeRD22 -0.62430 0.05636 -11.077 < 2e-16 ***
## FlowActionPeriodDuring:StationCodeI80 -0.51777 0.06421 -8.063 1.02e-15 ***
## FlowActionPeriodAfter:StationCodeI80 -0.50146 0.05684 -8.822 < 2e-16 ***
## FlowActionPeriodDuring:StationCodeLIS 0.42332 0.06385 6.630 3.90e-11 ***
## FlowActionPeriodAfter:StationCodeLIS -0.52929 0.05589 -9.470 < 2e-16 ***
## FlowActionPeriodDuring:StationCodeSTTD 1.12064 0.06775 16.541 < 2e-16 ***
## FlowActionPeriodAfter:StationCodeSTTD -0.38714 0.06322 -6.123 1.02e-09 ***
## FlowActionPeriodDuring:StationCodeLIB 0.26405 0.06912 3.820 0.000136 ***
## FlowActionPeriodAfter:StationCodeLIB -0.75378 0.06336 -11.897 < 2e-16 ***
## FlowActionPeriodDuring:StationCodeRVB 0.08567 0.06416 1.335 0.181922
## FlowActionPeriodAfter:StationCodeRVB -0.69367 0.05617 -12.350 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.772 3.974 27.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.864 Deviance explained = 86.6%
## -REML = 1561.7 Scale est. = 0.13479 n = 3478
appraise(m_gam_cat)
k.check(m_gam_cat)
## k' edf k-index p-value
## s(DOY) 4 3.77179 0.9457395 0
draw(m_gam_cat, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_cat, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_cat))
Box.test(residuals(m_gam_cat), lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: residuals(m_gam_cat)
## X-squared = 8609.3, df = 20, p-value < 2.2e-16
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_cat <- as.formula("Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_cat <- gamm(
f_gam_cat,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_cat_ar1 <- gamm(
f_gam_cat,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat_ar2 <- gamm(
f_gam_cat,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat_ar3 <- gamm(
f_gam_cat,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_cat$lme,
m_gamm_cat_ar1$lme,
m_gamm_cat_ar2$lme,
m_gamm_cat_ar3$lme
)
## Model df AIC BIC logLik Test L.Ratio
## m_gamm_cat$lme 1 60 3243.372 3611.616 -1561.686
## m_gamm_cat_ar1$lme 2 61 -2167.534 -1793.153 1144.767 1 vs 2 5412.906
## m_gamm_cat_ar2$lme 3 62 -2189.412 -1808.893 1156.706 2 vs 3 23.878
## m_gamm_cat_ar3$lme 4 63 -2188.415 -1801.759 1157.208 3 vs 4 1.004
## p-value
## m_gamm_cat$lme
## m_gamm_cat_ar1$lme <.0001
## m_gamm_cat_ar2$lme <.0001
## m_gamm_cat_ar3$lme 0.3164
It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
# Pull out the residuals and the GAM model
resid_cat_ar2 <- residuals(m_gamm_cat_ar2$lme, type = "normalized")
m_gamm_cat_ar2_gam <- m_gamm_cat_ar2$gam
summary(m_gamm_cat_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.092494 0.347301 26.180 < 2e-16
## Year_fct2016 0.074100 0.454595 0.163 0.870526
## Year_fct2017 0.351443 0.428540 0.820 0.412220
## Year_fct2018 0.220090 0.419794 0.524 0.600117
## Year_fct2019 -0.004065 0.421429 -0.010 0.992305
## FlowActionPeriodDuring 0.106229 0.094542 1.124 0.261253
## FlowActionPeriodAfter 0.205035 0.143274 1.431 0.152500
## StationCodeRD22 0.213335 0.427422 0.499 0.617727
## StationCodeI80 0.541690 0.439999 1.231 0.218364
## StationCodeLIS -0.597751 0.419022 -1.427 0.153805
## StationCodeSTTD -0.309144 0.429632 -0.720 0.471848
## StationCodeLIB -1.392380 0.421407 -3.304 0.000963
## StationCodeRVB -1.450049 0.417001 -3.477 0.000513
## Year_fct2016:FlowActionPeriodDuring -0.162149 0.090727 -1.787 0.073991
## Year_fct2017:FlowActionPeriodDuring 0.011836 0.089175 0.133 0.894417
## Year_fct2018:FlowActionPeriodDuring -0.171794 0.089080 -1.929 0.053872
## Year_fct2019:FlowActionPeriodDuring 0.054852 0.089154 0.615 0.538428
## Year_fct2016:FlowActionPeriodAfter -0.216376 0.127679 -1.695 0.090227
## Year_fct2017:FlowActionPeriodAfter -0.082755 0.125185 -0.661 0.508619
## Year_fct2018:FlowActionPeriodAfter -0.235315 0.125043 -1.882 0.059938
## Year_fct2019:FlowActionPeriodAfter 0.127736 0.124974 1.022 0.306806
## Year_fct2016:StationCodeRD22 -0.409097 0.576632 -0.709 0.478087
## Year_fct2017:StationCodeRD22 -0.510460 0.546951 -0.933 0.350740
## Year_fct2018:StationCodeRD22 -0.375917 0.536410 -0.701 0.483474
## Year_fct2019:StationCodeRD22 -0.125667 0.538556 -0.233 0.815511
## Year_fct2016:StationCodeI80 0.159128 0.585194 0.272 0.785697
## Year_fct2017:StationCodeI80 -0.229520 0.556103 -0.413 0.679830
## Year_fct2018:StationCodeI80 -0.577356 0.545768 -1.058 0.290186
## Year_fct2019:StationCodeI80 0.412105 0.547880 0.752 0.451994
## Year_fct2016:StationCodeLIS 0.766734 0.561958 1.364 0.172532
## Year_fct2017:StationCodeLIS 0.178449 0.544443 0.328 0.743110
## Year_fct2018:StationCodeLIS -0.130811 0.531634 -0.246 0.805655
## Year_fct2019:StationCodeLIS 0.384296 0.537308 0.715 0.474519
## Year_fct2016:StationCodeSTTD 0.494794 0.579473 0.854 0.393237
## Year_fct2017:StationCodeSTTD -0.666965 0.584959 -1.140 0.254287
## Year_fct2018:StationCodeSTTD -0.891467 0.556482 -1.602 0.109255
## Year_fct2019:StationCodeSTTD -1.068222 0.547691 -1.950 0.051209
## Year_fct2016:StationCodeLIB 1.032582 0.565695 1.825 0.068038
## Year_fct2017:StationCodeLIB -0.298959 0.549421 -0.544 0.586385
## Year_fct2018:StationCodeLIB -2.012752 0.597348 -3.369 0.000761
## Year_fct2019:StationCodeLIB -0.761233 0.605732 -1.257 0.208942
## Year_fct2016:StationCodeRVB 0.385100 0.575289 0.669 0.503284
## Year_fct2017:StationCodeRVB -0.578913 0.539279 -1.073 0.283125
## Year_fct2018:StationCodeRVB -0.447074 0.528583 -0.846 0.397726
## Year_fct2019:StationCodeRVB -0.854363 0.537540 -1.589 0.112064
## FlowActionPeriodDuring:StationCodeRD22 -0.067302 0.105406 -0.638 0.523192
## FlowActionPeriodAfter:StationCodeRD22 0.031609 0.148699 0.213 0.831675
## FlowActionPeriodDuring:StationCodeI80 -0.058480 0.105829 -0.553 0.580579
## FlowActionPeriodAfter:StationCodeI80 -0.002189 0.149048 -0.015 0.988284
## FlowActionPeriodDuring:StationCodeLIS -0.087621 0.105289 -0.832 0.405357
## FlowActionPeriodAfter:StationCodeLIS 0.013527 0.148521 0.091 0.927435
## FlowActionPeriodDuring:StationCodeSTTD 0.056894 0.105964 0.537 0.591360
## FlowActionPeriodAfter:StationCodeSTTD -0.493124 0.150394 -3.279 0.001053
## FlowActionPeriodDuring:StationCodeLIB -0.098209 0.106665 -0.921 0.357264
## FlowActionPeriodAfter:StationCodeLIB -0.224653 0.150868 -1.489 0.136561
## FlowActionPeriodDuring:StationCodeRVB -0.067906 0.105332 -0.645 0.519175
## FlowActionPeriodAfter:StationCodeRVB -0.174298 0.148565 -1.173 0.240792
##
## (Intercept) ***
## Year_fct2016
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring
## FlowActionPeriodAfter
## StationCodeRD22
## StationCodeI80
## StationCodeLIS
## StationCodeSTTD
## StationCodeLIB ***
## StationCodeRVB ***
## Year_fct2016:FlowActionPeriodDuring .
## Year_fct2017:FlowActionPeriodDuring
## Year_fct2018:FlowActionPeriodDuring .
## Year_fct2019:FlowActionPeriodDuring
## Year_fct2016:FlowActionPeriodAfter .
## Year_fct2017:FlowActionPeriodAfter
## Year_fct2018:FlowActionPeriodAfter .
## Year_fct2019:FlowActionPeriodAfter
## Year_fct2016:StationCodeRD22
## Year_fct2017:StationCodeRD22
## Year_fct2018:StationCodeRD22
## Year_fct2019:StationCodeRD22
## Year_fct2016:StationCodeI80
## Year_fct2017:StationCodeI80
## Year_fct2018:StationCodeI80
## Year_fct2019:StationCodeI80
## Year_fct2016:StationCodeLIS
## Year_fct2017:StationCodeLIS
## Year_fct2018:StationCodeLIS
## Year_fct2019:StationCodeLIS
## Year_fct2016:StationCodeSTTD
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD
## Year_fct2019:StationCodeSTTD .
## Year_fct2016:StationCodeLIB .
## Year_fct2017:StationCodeLIB
## Year_fct2018:StationCodeLIB ***
## Year_fct2019:StationCodeLIB
## Year_fct2016:StationCodeRVB
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## FlowActionPeriodDuring:StationCodeRD22
## FlowActionPeriodAfter:StationCodeRD22
## FlowActionPeriodDuring:StationCodeI80
## FlowActionPeriodAfter:StationCodeI80
## FlowActionPeriodDuring:StationCodeLIS
## FlowActionPeriodAfter:StationCodeLIS
## FlowActionPeriodDuring:StationCodeSTTD
## FlowActionPeriodAfter:StationCodeSTTD **
## FlowActionPeriodDuring:StationCodeLIB
## FlowActionPeriodAfter:StationCodeLIB
## FlowActionPeriodDuring:StationCodeRVB
## FlowActionPeriodAfter:StationCodeRVB
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.734 3.734 11.49 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.785
## Scale est. = 0.27538 n = 3478
appraise(m_gamm_cat_ar2_gam)
k.check(m_gamm_cat_ar2_gam)
## k' edf k-index p-value
## s(DOY) 4 3.734167 0.7698425 0
draw(m_gamm_cat_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_cat_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_cat_ar2)
Box.test(resid_cat_ar2, lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: resid_cat_ar2
## X-squared = 32.886, df = 20, p-value = 0.03473
The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_cat_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode)^2 + s(DOY,
## k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 0.345 0.8477
## FlowActionPeriod 2 1.070 0.3432
## StationCode 6 9.384 3.08e-10
## Year_fct:FlowActionPeriod 8 1.991 0.0438
## Year_fct:StationCode 24 2.732 1.16e-05
## FlowActionPeriod:StationCode 12 4.817 6.68e-08
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.734 3.734 11.49 <2e-16
All the interaction terms are significant in this model, so we’ll use this one for the model selection process.
We’ll try running the categorical model adding water temperature to see if it improves the predictive power of the model.
m_gam_cat_wt <- gam(
Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)
Lets look at the model summary and diagnostics:
summary(m_gam_cat_wt)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 +
## s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.4732371 0.4077027 23.236 < 2e-16
## Year_fct2016 0.1185304 0.3872590 0.306 0.759566
## Year_fct2017 0.1743035 0.3365746 0.518 0.604579
## Year_fct2018 0.3038120 0.2848691 1.066 0.286275
## Year_fct2019 0.0558398 0.2655069 0.210 0.833435
## FlowActionPeriodDuring -0.8448474 0.3953310 -2.137 0.032664
## FlowActionPeriodAfter 0.1409841 0.3375317 0.418 0.676199
## StationCodeRD22 -0.0797142 0.3149937 -0.253 0.800233
## StationCodeI80 1.0484912 0.3010805 3.482 0.000503
## StationCodeLIS -3.0959341 0.2998499 -10.325 < 2e-16
## StationCodeSTTD -2.1712782 0.3221696 -6.740 1.86e-11
## StationCodeLIB -2.4195633 0.3868708 -6.254 4.49e-10
## StationCodeRVB -2.5192484 0.3578309 -7.040 2.31e-12
## WaterTemp -0.0172525 0.0159383 -1.082 0.279127
## Year_fct2016:FlowActionPeriodDuring -0.1478412 0.0717935 -2.059 0.039545
## Year_fct2017:FlowActionPeriodDuring -0.1142243 0.0559969 -2.040 0.041444
## Year_fct2018:FlowActionPeriodDuring 0.0081453 0.0566778 0.144 0.885736
## Year_fct2019:FlowActionPeriodDuring -0.4834680 0.0537729 -8.991 < 2e-16
## Year_fct2016:FlowActionPeriodAfter -0.0463716 0.0939445 -0.494 0.621616
## Year_fct2017:FlowActionPeriodAfter -0.1164966 0.0985549 -1.182 0.237269
## Year_fct2018:FlowActionPeriodAfter -0.1779551 0.0818456 -2.174 0.029753
## Year_fct2019:FlowActionPeriodAfter 0.1434797 0.0868389 1.652 0.098575
## Year_fct2016:StationCodeRD22 -0.4030850 0.1061711 -3.797 0.000149
## Year_fct2017:StationCodeRD22 -0.6258445 0.0886093 -7.063 1.97e-12
## Year_fct2018:StationCodeRD22 -0.4032341 0.0861352 -4.681 2.96e-06
## Year_fct2019:StationCodeRD22 -0.2105029 0.0863540 -2.438 0.014833
## Year_fct2016:StationCodeI80 0.3764403 0.1078367 3.491 0.000488
## Year_fct2017:StationCodeI80 -0.2152108 0.0907540 -2.371 0.017778
## Year_fct2018:StationCodeI80 -0.5036823 0.0880693 -5.719 1.16e-08
## Year_fct2019:StationCodeI80 0.4707630 0.0885154 5.318 1.11e-07
## Year_fct2016:StationCodeLIS 0.3980072 0.1037768 3.835 0.000128
## Year_fct2017:StationCodeLIS 0.1331875 0.0890675 1.495 0.134915
## Year_fct2018:StationCodeLIS -0.0665137 0.0856235 -0.777 0.437321
## Year_fct2019:StationCodeLIS 0.5866469 0.0867510 6.762 1.59e-11
## Year_fct2016:StationCodeSTTD 0.4647981 0.1075972 4.320 1.61e-05
## Year_fct2017:StationCodeSTTD -0.1040733 0.0993597 -1.047 0.294971
## Year_fct2018:StationCodeSTTD -0.5941584 0.0894294 -6.644 3.54e-11
## Year_fct2019:StationCodeSTTD -0.8406972 0.0880556 -9.547 < 2e-16
## Year_fct2016:StationCodeLIB 1.0594550 0.1104642 9.591 < 2e-16
## Year_fct2017:StationCodeLIB -0.2288879 0.0930776 -2.459 0.013978
## Year_fct2018:StationCodeLIB -1.8694287 0.1011190 -18.487 < 2e-16
## Year_fct2019:StationCodeLIB -0.5687776 0.1049190 -5.421 6.34e-08
## Year_fct2016:StationCodeRVB 0.2273758 0.1104979 2.058 0.039691
## Year_fct2017:StationCodeRVB -0.4993343 0.0911536 -5.478 4.62e-08
## Year_fct2018:StationCodeRVB -0.3330771 0.0868161 -3.837 0.000127
## Year_fct2019:StationCodeRVB -0.5914447 0.0887079 -6.667 3.03e-11
## Year_fct2016:WaterTemp -0.0047186 0.0147895 -0.319 0.749709
## Year_fct2017:WaterTemp 0.0102785 0.0131868 0.779 0.435767
## Year_fct2018:WaterTemp -0.0074263 0.0113164 -0.656 0.511709
## Year_fct2019:WaterTemp -0.0009458 0.0102805 -0.092 0.926706
## FlowActionPeriodDuring:StationCodeRD22 -0.3367835 0.0692219 -4.865 1.19e-06
## FlowActionPeriodAfter:StationCodeRD22 -0.3952011 0.1035785 -3.815 0.000138
## FlowActionPeriodDuring:StationCodeI80 -0.5058395 0.0700402 -7.222 6.28e-13
## FlowActionPeriodAfter:StationCodeI80 -0.5403872 0.1029498 -5.249 1.62e-07
## FlowActionPeriodDuring:StationCodeLIS 0.6066616 0.0696598 8.709 < 2e-16
## FlowActionPeriodAfter:StationCodeLIS 0.1670770 0.0997364 1.675 0.093989
## FlowActionPeriodDuring:StationCodeSTTD 1.2258165 0.0734088 16.699 < 2e-16
## FlowActionPeriodAfter:StationCodeSTTD 0.0006864 0.1043389 0.007 0.994751
## FlowActionPeriodDuring:StationCodeLIB 0.4149955 0.0840477 4.938 8.29e-07
## FlowActionPeriodAfter:StationCodeLIB -0.4732155 0.1097713 -4.311 1.67e-05
## FlowActionPeriodDuring:StationCodeRVB 0.2264187 0.0775096 2.921 0.003510
## FlowActionPeriodAfter:StationCodeRVB -0.3932928 0.1064178 -3.696 0.000223
## FlowActionPeriodDuring:WaterTemp 0.0304368 0.0159083 1.913 0.055797
## FlowActionPeriodAfter:WaterTemp 0.0026541 0.0125480 0.212 0.832498
## StationCodeRD22:WaterTemp 0.0307775 0.0121747 2.528 0.011517
## StationCodeI80:WaterTemp -0.0098427 0.0115431 -0.853 0.393890
## StationCodeLIS:WaterTemp 0.1106623 0.0116009 9.539 < 2e-16
## StationCodeSTTD:WaterTemp 0.0615658 0.0127363 4.834 1.40e-06
## StationCodeLIB:WaterTemp 0.0474936 0.0159554 2.977 0.002935
## StationCodeRVB:WaterTemp 0.0499363 0.0143724 3.474 0.000518
##
## (Intercept) ***
## Year_fct2016
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring *
## FlowActionPeriodAfter
## StationCodeRD22
## StationCodeI80 ***
## StationCodeLIS ***
## StationCodeSTTD ***
## StationCodeLIB ***
## StationCodeRVB ***
## WaterTemp
## Year_fct2016:FlowActionPeriodDuring *
## Year_fct2017:FlowActionPeriodDuring *
## Year_fct2018:FlowActionPeriodDuring
## Year_fct2019:FlowActionPeriodDuring ***
## Year_fct2016:FlowActionPeriodAfter
## Year_fct2017:FlowActionPeriodAfter
## Year_fct2018:FlowActionPeriodAfter *
## Year_fct2019:FlowActionPeriodAfter .
## Year_fct2016:StationCodeRD22 ***
## Year_fct2017:StationCodeRD22 ***
## Year_fct2018:StationCodeRD22 ***
## Year_fct2019:StationCodeRD22 *
## Year_fct2016:StationCodeI80 ***
## Year_fct2017:StationCodeI80 *
## Year_fct2018:StationCodeI80 ***
## Year_fct2019:StationCodeI80 ***
## Year_fct2016:StationCodeLIS ***
## Year_fct2017:StationCodeLIS
## Year_fct2018:StationCodeLIS
## Year_fct2019:StationCodeLIS ***
## Year_fct2016:StationCodeSTTD ***
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD ***
## Year_fct2019:StationCodeSTTD ***
## Year_fct2016:StationCodeLIB ***
## Year_fct2017:StationCodeLIB *
## Year_fct2018:StationCodeLIB ***
## Year_fct2019:StationCodeLIB ***
## Year_fct2016:StationCodeRVB *
## Year_fct2017:StationCodeRVB ***
## Year_fct2018:StationCodeRVB ***
## Year_fct2019:StationCodeRVB ***
## Year_fct2016:WaterTemp
## Year_fct2017:WaterTemp
## Year_fct2018:WaterTemp
## Year_fct2019:WaterTemp
## FlowActionPeriodDuring:StationCodeRD22 ***
## FlowActionPeriodAfter:StationCodeRD22 ***
## FlowActionPeriodDuring:StationCodeI80 ***
## FlowActionPeriodAfter:StationCodeI80 ***
## FlowActionPeriodDuring:StationCodeLIS ***
## FlowActionPeriodAfter:StationCodeLIS .
## FlowActionPeriodDuring:StationCodeSTTD ***
## FlowActionPeriodAfter:StationCodeSTTD
## FlowActionPeriodDuring:StationCodeLIB ***
## FlowActionPeriodAfter:StationCodeLIB ***
## FlowActionPeriodDuring:StationCodeRVB **
## FlowActionPeriodAfter:StationCodeRVB ***
## FlowActionPeriodDuring:WaterTemp .
## FlowActionPeriodAfter:WaterTemp
## StationCodeRD22:WaterTemp *
## StationCodeI80:WaterTemp
## StationCodeLIS:WaterTemp ***
## StationCodeSTTD:WaterTemp ***
## StationCodeLIB:WaterTemp **
## StationCodeRVB:WaterTemp ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.731 3.966 23.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.87 Deviance explained = 87.3%
## -REML = 1517.2 Scale est. = 0.12822 n = 3478
appraise(m_gam_cat_wt)
k.check(m_gam_cat_wt)
## k' edf k-index p-value
## s(DOY) 4 3.730833 0.9463198 0
draw(m_gam_cat_wt, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_cat_wt, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_cat_wt))
Box.test(residuals(m_gam_cat_wt), lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: residuals(m_gam_cat_wt)
## X-squared = 8031.8, df = 20, p-value < 2.2e-16
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_cat_wt <- as.formula("Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_cat_wt <- gamm(
f_gam_cat_wt,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_cat_wt_ar1 <- gamm(
f_gam_cat_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat_wt_ar2 <- gamm(
f_gam_cat_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_cat_wt_ar3 <- gamm(
f_gam_cat_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_cat_wt$lme,
m_gamm_cat_wt_ar1$lme,
m_gamm_cat_wt_ar2$lme,
m_gamm_cat_wt_ar3$lme
)
## Model df AIC BIC logLik Test L.Ratio
## m_gamm_cat_wt$lme 1 73 3180.470 3628.222 -1517.235
## m_gamm_cat_wt_ar1$lme 2 74 -2135.571 -1681.686 1141.786 1 vs 2 5318.041
## m_gamm_cat_wt_ar2$lme 3 75 -2160.486 -1700.467 1155.243 2 vs 3 26.915
## m_gamm_cat_wt_ar3$lme 4 76 -2159.702 -1693.550 1155.851 3 vs 4 1.216
## p-value
## m_gamm_cat_wt$lme
## m_gamm_cat_wt_ar1$lme <.0001
## m_gamm_cat_wt_ar2$lme <.0001
## m_gamm_cat_wt_ar3$lme 0.2701
It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
# Pull out the residuals and the GAM model
resid_cat_wt_ar2 <- residuals(m_gamm_cat_wt_ar2$lme, type = "normalized")
m_gamm_cat_wt_ar2_gam <- m_gamm_cat_wt_ar2$gam
summary(m_gamm_cat_wt_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 +
## s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.1445090 0.6901525 13.250 < 2e-16
## Year_fct2016 1.3780976 0.6921213 1.991 0.04655
## Year_fct2017 0.0579152 0.6440942 0.090 0.92836
## Year_fct2018 -0.3370015 0.6545952 -0.515 0.60671
## Year_fct2019 -0.3900285 0.6187327 -0.630 0.52850
## FlowActionPeriodDuring -0.4663337 0.4240927 -1.100 0.27158
## FlowActionPeriodAfter -0.3836823 0.3632026 -1.056 0.29087
## StationCodeRD22 -1.1868159 0.6842598 -1.734 0.08293
## StationCodeI80 0.9164038 0.6678780 1.372 0.17012
## StationCodeLIS -1.6220417 0.6764787 -2.398 0.01655
## StationCodeSTTD -0.5586425 0.6987943 -0.799 0.42409
## StationCodeLIB -0.2983707 0.7674643 -0.389 0.69747
## StationCodeRVB -1.9098221 0.8359931 -2.284 0.02240
## WaterTemp 0.0009966 0.0243656 0.041 0.96738
## Year_fct2016:FlowActionPeriodDuring -0.1572496 0.0926170 -1.698 0.08963
## Year_fct2017:FlowActionPeriodDuring -0.0258551 0.0919774 -0.281 0.77865
## Year_fct2018:FlowActionPeriodDuring -0.1036021 0.0974311 -1.063 0.28770
## Year_fct2019:FlowActionPeriodDuring 0.0270928 0.0904924 0.299 0.76466
## Year_fct2016:FlowActionPeriodAfter -0.1711207 0.1383100 -1.237 0.21609
## Year_fct2017:FlowActionPeriodAfter -0.0309569 0.1287914 -0.240 0.81006
## Year_fct2018:FlowActionPeriodAfter -0.1354965 0.1343868 -1.008 0.31340
## Year_fct2019:FlowActionPeriodAfter 0.1469567 0.1263950 1.163 0.24504
## Year_fct2016:StationCodeRD22 -0.6182502 0.6279173 -0.985 0.32489
## Year_fct2017:StationCodeRD22 -0.5837985 0.5953210 -0.981 0.32684
## Year_fct2018:StationCodeRD22 -0.3239681 0.5851972 -0.554 0.57989
## Year_fct2019:StationCodeRD22 -0.1547261 0.5869605 -0.264 0.79210
## Year_fct2016:StationCodeI80 0.2258097 0.6387648 0.354 0.72373
## Year_fct2017:StationCodeI80 -0.1772145 0.6064542 -0.292 0.77014
## Year_fct2018:StationCodeI80 -0.5434617 0.5963020 -0.911 0.36216
## Year_fct2019:StationCodeI80 0.4333058 0.5981682 0.724 0.46888
## Year_fct2016:StationCodeLIS 0.5871339 0.6122516 0.959 0.33764
## Year_fct2017:StationCodeLIS 0.1547944 0.5925727 0.261 0.79394
## Year_fct2018:StationCodeLIS -0.1245806 0.5799582 -0.215 0.82993
## Year_fct2019:StationCodeLIS 0.3524021 0.5857178 0.602 0.54744
## Year_fct2016:StationCodeSTTD 0.4354888 0.6320048 0.689 0.49083
## Year_fct2017:StationCodeSTTD -0.6528634 0.6368032 -1.025 0.30533
## Year_fct2018:StationCodeSTTD -0.8898961 0.6067842 -1.467 0.14258
## Year_fct2019:StationCodeSTTD -1.0781623 0.5972945 -1.805 0.07115
## Year_fct2016:StationCodeLIB 0.9323841 0.6169532 1.511 0.13081
## Year_fct2017:StationCodeLIB -0.2392731 0.5986330 -0.400 0.68940
## Year_fct2018:StationCodeLIB -2.0226118 0.6499257 -3.112 0.00187
## Year_fct2019:StationCodeLIB -0.7318547 0.6572217 -1.114 0.26555
## Year_fct2016:StationCodeRVB 0.2238684 0.6270456 0.357 0.72110
## Year_fct2017:StationCodeRVB -0.5138895 0.5877741 -0.874 0.38202
## Year_fct2018:StationCodeRVB -0.3928585 0.5778469 -0.680 0.49664
## Year_fct2019:StationCodeRVB -0.7921245 0.5874493 -1.348 0.17762
## Year_fct2016:WaterTemp -0.0524611 0.0207223 -2.532 0.01140
## Year_fct2017:WaterTemp 0.0126743 0.0191452 0.662 0.50801
## Year_fct2018:WaterTemp 0.0245102 0.0205942 1.190 0.23407
## Year_fct2019:WaterTemp 0.0193582 0.0181899 1.064 0.28730
## FlowActionPeriodDuring:StationCodeRD22 -0.0578678 0.1049823 -0.551 0.58152
## FlowActionPeriodAfter:StationCodeRD22 0.0980816 0.1507379 0.651 0.51530
## FlowActionPeriodDuring:StationCodeI80 -0.0674505 0.1057957 -0.638 0.52381
## FlowActionPeriodAfter:StationCodeI80 -0.0037150 0.1519709 -0.024 0.98050
## FlowActionPeriodDuring:StationCodeLIS -0.0502238 0.1055792 -0.476 0.63432
## FlowActionPeriodAfter:StationCodeLIS 0.1157942 0.1515187 0.764 0.44479
## FlowActionPeriodDuring:StationCodeSTTD 0.0666638 0.1064582 0.626 0.53123
## FlowActionPeriodAfter:StationCodeSTTD -0.4517271 0.1532547 -2.948 0.00322
## FlowActionPeriodDuring:StationCodeLIB -0.0730926 0.1151386 -0.635 0.52559
## FlowActionPeriodAfter:StationCodeLIB -0.2143517 0.1567062 -1.368 0.17145
## FlowActionPeriodDuring:StationCodeRVB -0.0110737 0.1127300 -0.098 0.92175
## FlowActionPeriodAfter:StationCodeRVB -0.0964396 0.1546565 -0.624 0.53295
## FlowActionPeriodDuring:WaterTemp 0.0231542 0.0168338 1.375 0.16908
## FlowActionPeriodAfter:WaterTemp 0.0210524 0.0143386 1.468 0.14213
## StationCodeRD22:WaterTemp 0.0653478 0.0222444 2.938 0.00333
## StationCodeI80:WaterTemp -0.0191793 0.0207291 -0.925 0.35491
## StationCodeLIS:WaterTemp 0.0464759 0.0218848 2.124 0.03377
## StationCodeSTTD:WaterTemp 0.0103409 0.0232382 0.445 0.65635
## StationCodeLIB:WaterTemp -0.0551245 0.0276723 -1.992 0.04645
## StationCodeRVB:WaterTemp 0.0192565 0.0319521 0.603 0.54677
##
## (Intercept) ***
## Year_fct2016 *
## Year_fct2017
## Year_fct2018
## Year_fct2019
## FlowActionPeriodDuring
## FlowActionPeriodAfter
## StationCodeRD22 .
## StationCodeI80
## StationCodeLIS *
## StationCodeSTTD
## StationCodeLIB
## StationCodeRVB *
## WaterTemp
## Year_fct2016:FlowActionPeriodDuring .
## Year_fct2017:FlowActionPeriodDuring
## Year_fct2018:FlowActionPeriodDuring
## Year_fct2019:FlowActionPeriodDuring
## Year_fct2016:FlowActionPeriodAfter
## Year_fct2017:FlowActionPeriodAfter
## Year_fct2018:FlowActionPeriodAfter
## Year_fct2019:FlowActionPeriodAfter
## Year_fct2016:StationCodeRD22
## Year_fct2017:StationCodeRD22
## Year_fct2018:StationCodeRD22
## Year_fct2019:StationCodeRD22
## Year_fct2016:StationCodeI80
## Year_fct2017:StationCodeI80
## Year_fct2018:StationCodeI80
## Year_fct2019:StationCodeI80
## Year_fct2016:StationCodeLIS
## Year_fct2017:StationCodeLIS
## Year_fct2018:StationCodeLIS
## Year_fct2019:StationCodeLIS
## Year_fct2016:StationCodeSTTD
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD
## Year_fct2019:StationCodeSTTD .
## Year_fct2016:StationCodeLIB
## Year_fct2017:StationCodeLIB
## Year_fct2018:StationCodeLIB **
## Year_fct2019:StationCodeLIB
## Year_fct2016:StationCodeRVB
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## Year_fct2016:WaterTemp *
## Year_fct2017:WaterTemp
## Year_fct2018:WaterTemp
## Year_fct2019:WaterTemp
## FlowActionPeriodDuring:StationCodeRD22
## FlowActionPeriodAfter:StationCodeRD22
## FlowActionPeriodDuring:StationCodeI80
## FlowActionPeriodAfter:StationCodeI80
## FlowActionPeriodDuring:StationCodeLIS
## FlowActionPeriodAfter:StationCodeLIS
## FlowActionPeriodDuring:StationCodeSTTD
## FlowActionPeriodAfter:StationCodeSTTD **
## FlowActionPeriodDuring:StationCodeLIB
## FlowActionPeriodAfter:StationCodeLIB
## FlowActionPeriodDuring:StationCodeRVB
## FlowActionPeriodAfter:StationCodeRVB
## FlowActionPeriodDuring:WaterTemp
## FlowActionPeriodAfter:WaterTemp
## StationCodeRD22:WaterTemp **
## StationCodeI80:WaterTemp
## StationCodeLIS:WaterTemp *
## StationCodeSTTD:WaterTemp
## StationCodeLIB:WaterTemp *
## StationCodeRVB:WaterTemp
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.669 3.669 8.849 5.92e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.774
## Scale est. = 0.30258 n = 3478
appraise(m_gamm_cat_wt_ar2_gam)
k.check(m_gamm_cat_wt_ar2_gam)
## k' edf k-index p-value
## s(DOY) 4 3.669285 0.7410218 0
draw(m_gamm_cat_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_cat_wt_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_cat_wt_ar2)
Box.test(resid_cat_wt_ar2, lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: resid_cat_wt_ar2
## X-squared = 30.859, df = 20, p-value = 0.05708
The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_cat_wt_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + FlowActionPeriod + StationCode + WaterTemp)^2 +
## s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 2.412 0.047021
## FlowActionPeriod 2 0.740 0.477314
## StationCode 6 5.258 2.11e-05
## WaterTemp 1 0.002 0.967377
## Year_fct:FlowActionPeriod 8 1.159 0.320505
## Year_fct:StationCode 24 2.285 0.000357
## Year_fct:WaterTemp 4 4.446 0.001386
## FlowActionPeriod:StationCode 12 4.941 3.63e-08
## FlowActionPeriod:WaterTemp 2 1.414 0.243334
## StationCode:WaterTemp 6 6.427 9.53e-07
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.669 3.669 8.849 5.92e-06
The Year:StationCode, Year:WaterTemp, FlowActionPeriod:StationCode, and StationCode:WaterTemp interactions were significant among the interaction terms of the model. Let’s remove the Year:FlowActionPeriod and FlowActionPeriod:WaterTemp interaction terms and check the model summary and diagnostics.
m_gamm_cat_wt_ar2b <- gamm(
Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode * FlowActionPeriod + s(DOY, k = 5),
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
# Pull out the residuals and the GAM model
resid_cat_wt_ar2b <- residuals(m_gamm_cat_wt_ar2b$lme, type = "normalized")
m_gamm_cat_wt_ar2b_gam <- m_gamm_cat_wt_ar2b$gam
summary(m_gamm_cat_wt_ar2b_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode *
## FlowActionPeriod + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.782805 0.664490 13.217 < 2e-16
## Year_fct2016 1.341143 0.707272 1.896 0.05802
## Year_fct2017 -0.006022 0.652080 -0.009 0.99263
## Year_fct2018 -0.528205 0.657513 -0.803 0.42184
## Year_fct2019 -0.272579 0.628678 -0.434 0.66462
## StationCodeRD22 -1.163183 0.707102 -1.645 0.10006
## StationCodeI80 0.951126 0.691075 1.376 0.16882
## StationCodeLIS -1.559134 0.698191 -2.233 0.02561
## StationCodeSTTD -0.493857 0.720065 -0.686 0.49285
## StationCodeLIB -0.085665 0.784497 -0.109 0.91305
## StationCodeRVB -1.782352 0.855977 -2.082 0.03739
## WaterTemp 0.016614 0.022904 0.725 0.46828
## FlowActionPeriodDuring 0.052204 0.074987 0.696 0.48637
## FlowActionPeriodAfter 0.099477 0.110745 0.898 0.36912
## Year_fct2016:StationCodeRD22 -0.613100 0.667140 -0.919 0.35816
## Year_fct2017:StationCodeRD22 -0.578281 0.633452 -0.913 0.36136
## Year_fct2018:StationCodeRD22 -0.315304 0.622783 -0.506 0.61269
## Year_fct2019:StationCodeRD22 -0.149630 0.624636 -0.240 0.81070
## Year_fct2016:StationCodeI80 0.223843 0.678341 0.330 0.74143
## Year_fct2017:StationCodeI80 -0.181161 0.644952 -0.281 0.77881
## Year_fct2018:StationCodeI80 -0.547381 0.634232 -0.863 0.38817
## Year_fct2019:StationCodeI80 0.423004 0.636191 0.665 0.50616
## Year_fct2016:StationCodeLIS 0.619479 0.650406 0.952 0.34094
## Year_fct2017:StationCodeLIS 0.170080 0.630490 0.270 0.78736
## Year_fct2018:StationCodeLIS -0.116806 0.617195 -0.189 0.84991
## Year_fct2019:StationCodeLIS 0.349457 0.623223 0.561 0.57502
## Year_fct2016:StationCodeSTTD 0.445998 0.671203 0.664 0.50643
## Year_fct2017:StationCodeSTTD -0.666675 0.676269 -0.986 0.32429
## Year_fct2018:StationCodeSTTD -0.889459 0.645302 -1.378 0.16818
## Year_fct2019:StationCodeSTTD -1.078389 0.635280 -1.698 0.08969
## Year_fct2016:StationCodeLIB 0.940927 0.655358 1.436 0.15117
## Year_fct2017:StationCodeLIB -0.223119 0.636768 -0.350 0.72607
## Year_fct2018:StationCodeLIB -2.019727 0.689265 -2.930 0.00341
## Year_fct2019:StationCodeLIB -0.746492 0.696748 -1.071 0.28407
## Year_fct2016:StationCodeRVB 0.239251 0.665883 0.359 0.71939
## Year_fct2017:StationCodeRVB -0.506117 0.625396 -0.809 0.41841
## Year_fct2018:StationCodeRVB -0.384654 0.614812 -0.626 0.53159
## Year_fct2019:StationCodeRVB -0.807469 0.624907 -1.292 0.19640
## Year_fct2016:WaterTemp -0.056611 0.020507 -2.761 0.00580
## Year_fct2017:WaterTemp 0.013897 0.018782 0.740 0.45940
## Year_fct2018:WaterTemp 0.029403 0.020156 1.459 0.14472
## Year_fct2019:WaterTemp 0.016339 0.017903 0.913 0.36148
## StationCodeRD22:WaterTemp 0.064117 0.022486 2.851 0.00438
## StationCodeI80:WaterTemp -0.020588 0.020915 -0.984 0.32500
## StationCodeLIS:WaterTemp 0.043491 0.022089 1.969 0.04904
## StationCodeSTTD:WaterTemp 0.008117 0.023400 0.347 0.72870
## StationCodeLIB:WaterTemp -0.063750 0.027815 -2.292 0.02197
## StationCodeRVB:WaterTemp 0.015016 0.032392 0.464 0.64298
## StationCodeRD22:FlowActionPeriodDuring -0.054000 0.105338 -0.513 0.60824
## StationCodeI80:FlowActionPeriodDuring -0.062165 0.106067 -0.586 0.55785
## StationCodeLIS:FlowActionPeriodDuring -0.061674 0.105451 -0.585 0.55868
## StationCodeSTTD:FlowActionPeriodDuring 0.042061 0.105724 0.398 0.69078
## StationCodeLIB:FlowActionPeriodDuring -0.133669 0.107077 -1.248 0.21199
## StationCodeRVB:FlowActionPeriodDuring -0.063906 0.105456 -0.606 0.54456
## StationCodeRD22:FlowActionPeriodAfter 0.101187 0.150452 0.673 0.50128
## StationCodeI80:FlowActionPeriodAfter 0.005898 0.150783 0.039 0.96880
## StationCodeLIS:FlowActionPeriodAfter 0.107314 0.151165 0.710 0.47781
## StationCodeSTTD:FlowActionPeriodAfter -0.476380 0.152352 -3.127 0.00178
## StationCodeLIB:FlowActionPeriodAfter -0.265101 0.153142 -1.731 0.08353
## StationCodeRVB:FlowActionPeriodAfter -0.142648 0.151288 -0.943 0.34580
##
## (Intercept) ***
## Year_fct2016 .
## Year_fct2017
## Year_fct2018
## Year_fct2019
## StationCodeRD22
## StationCodeI80
## StationCodeLIS *
## StationCodeSTTD
## StationCodeLIB
## StationCodeRVB *
## WaterTemp
## FlowActionPeriodDuring
## FlowActionPeriodAfter
## Year_fct2016:StationCodeRD22
## Year_fct2017:StationCodeRD22
## Year_fct2018:StationCodeRD22
## Year_fct2019:StationCodeRD22
## Year_fct2016:StationCodeI80
## Year_fct2017:StationCodeI80
## Year_fct2018:StationCodeI80
## Year_fct2019:StationCodeI80
## Year_fct2016:StationCodeLIS
## Year_fct2017:StationCodeLIS
## Year_fct2018:StationCodeLIS
## Year_fct2019:StationCodeLIS
## Year_fct2016:StationCodeSTTD
## Year_fct2017:StationCodeSTTD
## Year_fct2018:StationCodeSTTD
## Year_fct2019:StationCodeSTTD .
## Year_fct2016:StationCodeLIB
## Year_fct2017:StationCodeLIB
## Year_fct2018:StationCodeLIB **
## Year_fct2019:StationCodeLIB
## Year_fct2016:StationCodeRVB
## Year_fct2017:StationCodeRVB
## Year_fct2018:StationCodeRVB
## Year_fct2019:StationCodeRVB
## Year_fct2016:WaterTemp **
## Year_fct2017:WaterTemp
## Year_fct2018:WaterTemp
## Year_fct2019:WaterTemp
## StationCodeRD22:WaterTemp **
## StationCodeI80:WaterTemp
## StationCodeLIS:WaterTemp *
## StationCodeSTTD:WaterTemp
## StationCodeLIB:WaterTemp *
## StationCodeRVB:WaterTemp
## StationCodeRD22:FlowActionPeriodDuring
## StationCodeI80:FlowActionPeriodDuring
## StationCodeLIS:FlowActionPeriodDuring
## StationCodeSTTD:FlowActionPeriodDuring
## StationCodeLIB:FlowActionPeriodDuring
## StationCodeRVB:FlowActionPeriodDuring
## StationCodeRD22:FlowActionPeriodAfter
## StationCodeI80:FlowActionPeriodAfter
## StationCodeLIS:FlowActionPeriodAfter
## StationCodeSTTD:FlowActionPeriodAfter **
## StationCodeLIB:FlowActionPeriodAfter .
## StationCodeRVB:FlowActionPeriodAfter
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.656 3.656 7.364 1.34e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.762
## Scale est. = 0.32529 n = 3478
appraise(m_gamm_cat_wt_ar2b_gam)
k.check(m_gamm_cat_wt_ar2b_gam)
## k' edf k-index p-value
## s(DOY) 4 3.655508 0.7254327 0
draw(m_gamm_cat_wt_ar2b_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_cat_wt_ar2b_gam, pages = 1, all.terms = TRUE)
acf(resid_cat_wt_ar2b)
Box.test(resid_cat_wt_ar2b, lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: resid_cat_wt_ar2b
## X-squared = 31.046, df = 20, p-value = 0.05459
anova(m_gamm_cat_wt_ar2b_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + StationCode + WaterTemp)^2 + StationCode *
## FlowActionPeriod + s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 2.316 0.055024
## StationCode 6 4.867 5.84e-05
## WaterTemp 1 0.526 0.468278
## FlowActionPeriod 2 0.415 0.660676
## Year_fct:StationCode 24 2.036 0.002086
## Year_fct:WaterTemp 4 5.222 0.000342
## StationCode:WaterTemp 6 6.644 5.32e-07
## StationCode:FlowActionPeriod 12 5.012 2.55e-08
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.656 3.656 7.364 1.34e-05
This model still looks pretty good. We’ll use this one for the model selection process.
df_chla_wt_q %>%
ggplot(aes(x = Flow, y = Chla)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(vars(StationCode), scales = "free") +
theme_bw()
Let’s break these scatterplots apart by year to see how these patterns vary annually.
df_chla_wt_q %>%
ggplot(aes(x = Flow, y = Chla)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(
vars(StationCode, Year),
ncol = 5,
scales = "free",
labeller = labeller(.multi_line = FALSE)
) +
theme_bw()
We’ll try running a GAM including all two-way interactions between Year, Daily Average Flow, and Station, and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_q <- gam(
Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)
Lets look at the model summary and diagnostics:
summary(m_gam_q)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.322e+00 6.781e-02 137.468 < 2e-16 ***
## Year_fct2016 -6.414e-02 7.625e-02 -0.841 0.400279
## Year_fct2017 1.818e-01 7.628e-02 2.383 0.017221 *
## Year_fct2018 1.125e-01 7.096e-02 1.585 0.113065
## Year_fct2019 1.372e-02 7.116e-02 0.193 0.847183
## Flow -1.275e-03 1.059e-04 -12.038 < 2e-16 ***
## StationCodeRD22 3.344e-01 7.804e-02 4.285 1.88e-05 ***
## StationCodeI80 5.559e-01 7.984e-02 6.963 3.98e-12 ***
## StationCodeLIS -6.816e-01 7.585e-02 -8.987 < 2e-16 ***
## StationCodeSTTD -6.447e-01 7.727e-02 -8.344 < 2e-16 ***
## StationCodeLIB -1.640e+00 1.260e-01 -13.019 < 2e-16 ***
## StationCodeRVB -1.908e+00 1.136e-01 -16.801 < 2e-16 ***
## Year_fct2016:Flow -8.355e-05 3.159e-05 -2.645 0.008212 **
## Year_fct2017:Flow -7.085e-05 2.628e-05 -2.696 0.007046 **
## Year_fct2018:Flow -2.494e-05 2.636e-05 -0.946 0.344219
## Year_fct2019:Flow -5.410e-05 2.800e-05 -1.932 0.053386 .
## Year_fct2016:StationCodeRD22 -4.235e-01 9.371e-02 -4.519 6.41e-06 ***
## Year_fct2017:StationCodeRD22 -5.668e-01 9.263e-02 -6.119 1.05e-09 ***
## Year_fct2018:StationCodeRD22 -4.231e-01 8.681e-02 -4.874 1.14e-06 ***
## Year_fct2019:StationCodeRD22 -1.754e-01 8.709e-02 -2.014 0.044067 *
## Year_fct2016:StationCodeI80 2.386e-01 9.525e-02 2.505 0.012300 *
## Year_fct2017:StationCodeI80 -2.309e-01 9.413e-02 -2.453 0.014198 *
## Year_fct2018:StationCodeI80 -4.800e-01 8.848e-02 -5.425 6.21e-08 ***
## Year_fct2019:StationCodeI80 5.216e-01 8.876e-02 5.876 4.60e-09 ***
## Year_fct2016:StationCodeLIS 5.173e-01 9.177e-02 5.637 1.87e-08 ***
## Year_fct2017:StationCodeLIS 3.056e-01 9.172e-02 3.332 0.000870 ***
## Year_fct2018:StationCodeLIS -1.288e-01 8.627e-02 -1.493 0.135495
## Year_fct2019:StationCodeLIS 5.536e-01 8.765e-02 6.316 3.02e-10 ***
## Year_fct2016:StationCodeSTTD 3.119e-01 9.445e-02 3.302 0.000971 ***
## Year_fct2017:StationCodeSTTD -1.372e-01 9.943e-02 -1.380 0.167655
## Year_fct2018:StationCodeSTTD -7.586e-01 9.080e-02 -8.354 < 2e-16 ***
## Year_fct2019:StationCodeSTTD -1.145e+00 8.934e-02 -12.815 < 2e-16 ***
## Year_fct2016:StationCodeLIB 8.329e-01 1.058e-01 7.871 4.67e-15 ***
## Year_fct2017:StationCodeLIB -1.933e-01 1.063e-01 -1.818 0.069088 .
## Year_fct2018:StationCodeLIB -1.979e+00 1.223e-01 -16.186 < 2e-16 ***
## Year_fct2019:StationCodeLIB -7.144e-01 1.241e-01 -5.757 9.31e-09 ***
## Year_fct2016:StationCodeRVB 6.335e-01 2.347e-01 2.699 0.006986 **
## Year_fct2017:StationCodeRVB -2.198e-02 1.987e-01 -0.111 0.911935
## Year_fct2018:StationCodeRVB -4.140e-01 1.626e-01 -2.546 0.010955 *
## Year_fct2019:StationCodeRVB -4.845e-01 1.910e-01 -2.537 0.011238 *
## Flow:StationCodeRD22 -1.057e-04 1.375e-04 -0.768 0.442446
## Flow:StationCodeI80 -6.752e-04 1.387e-04 -4.867 1.18e-06 ***
## Flow:StationCodeLIS 1.541e-03 1.391e-04 11.079 < 2e-16 ***
## Flow:StationCodeSTTD 3.211e-03 1.406e-04 22.832 < 2e-16 ***
## Flow:StationCodeLIB 1.275e-03 1.220e-04 10.452 < 2e-16 ***
## Flow:StationCodeRVB 1.327e-03 1.039e-04 12.773 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.757 3.966 81.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.851 Deviance explained = 85.3%
## -REML = 1770.4 Scale est. = 0.14725 n = 3478
appraise(m_gam_q)
k.check(m_gam_q)
## k' edf k-index p-value
## s(DOY) 4 3.756521 0.9228792 0
draw(m_gam_q, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_q, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_q))
Box.test(residuals(m_gam_q), lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: residuals(m_gam_q)
## X-squared = 10096, df = 20, p-value < 2.2e-16
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_q <- as.formula("Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q <- gamm(
f_gam_q,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_q_ar1 <- gamm(
f_gam_q,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_ar2 <- gamm(
f_gam_q,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_ar3 <- gamm(
f_gam_q,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_q$lme,
m_gamm_q_ar1$lme,
m_gamm_q_ar2$lme,
m_gamm_q_ar3$lme
)
## Model df AIC BIC logLik Test L.Ratio p-value
## m_gamm_q$lme 1 49 3638.827 3939.716 -1770.413
## m_gamm_q_ar1$lme 2 50 -2092.217 -1785.187 1096.109 1 vs 2 5733.044 <.0001
## m_gamm_q_ar2$lme 3 51 -2111.366 -1798.195 1106.683 2 vs 3 21.149 <.0001
## m_gamm_q_ar3$lme 4 52 -2110.010 -1790.698 1107.005 3 vs 4 0.644 0.4224
It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
# Pull out the residuals and the GAM model
resid_q_ar2 <- residuals(m_gamm_q_ar2$lme, type = "normalized")
m_gamm_q_ar2_gam <- m_gamm_q_ar2$gam
summary(m_gamm_q_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.178e+00 2.921e-01 31.417 < 2e-16 ***
## Year_fct2016 -3.270e-02 3.761e-01 -0.087 0.930719
## Year_fct2017 3.329e-01 3.587e-01 0.928 0.353319
## Year_fct2018 1.306e-01 3.495e-01 0.374 0.708602
## Year_fct2019 8.786e-02 3.505e-01 0.251 0.802084
## Flow -1.821e-04 1.772e-04 -1.028 0.304126
## StationCodeRD22 3.170e-01 3.570e-01 0.888 0.374697
## StationCodeI80 7.103e-01 3.671e-01 1.935 0.053071 .
## StationCodeLIS -5.832e-01 3.496e-01 -1.668 0.095356 .
## StationCodeSTTD -4.766e-01 3.576e-01 -1.333 0.182705
## StationCodeLIB -1.228e+00 3.569e-01 -3.440 0.000588 ***
## StationCodeRVB -1.517e+00 3.519e-01 -4.312 1.66e-05 ***
## Year_fct2016:Flow -5.265e-06 1.924e-05 -0.274 0.784344
## Year_fct2017:Flow -4.817e-06 1.964e-05 -0.245 0.806252
## Year_fct2018:Flow -3.172e-06 2.102e-05 -0.151 0.880067
## Year_fct2019:Flow -2.436e-05 1.961e-05 -1.242 0.214151
## Year_fct2016:StationCodeRD22 -4.412e-01 4.823e-01 -0.915 0.360450
## Year_fct2017:StationCodeRD22 -5.891e-01 4.591e-01 -1.283 0.199529
## Year_fct2018:StationCodeRD22 -4.262e-01 4.486e-01 -0.950 0.342068
## Year_fct2019:StationCodeRD22 -1.675e-01 4.503e-01 -0.372 0.709983
## Year_fct2016:StationCodeI80 1.284e-01 4.896e-01 0.262 0.793087
## Year_fct2017:StationCodeI80 -3.791e-01 4.669e-01 -0.812 0.416906
## Year_fct2018:StationCodeI80 -6.219e-01 4.564e-01 -1.363 0.173092
## Year_fct2019:StationCodeI80 3.949e-01 4.582e-01 0.862 0.388749
## Year_fct2016:StationCodeLIS 6.885e-01 4.700e-01 1.465 0.143014
## Year_fct2017:StationCodeLIS 1.647e-01 4.566e-01 0.361 0.718331
## Year_fct2018:StationCodeLIS -1.577e-01 4.445e-01 -0.355 0.722725
## Year_fct2019:StationCodeLIS 4.135e-01 4.496e-01 0.920 0.357786
## Year_fct2016:StationCodeSTTD 3.444e-01 4.849e-01 0.710 0.477618
## Year_fct2017:StationCodeSTTD -5.748e-01 4.914e-01 -1.170 0.242151
## Year_fct2018:StationCodeSTTD -8.724e-01 4.660e-01 -1.872 0.061252 .
## Year_fct2019:StationCodeSTTD -1.196e+00 4.584e-01 -2.609 0.009123 **
## Year_fct2016:StationCodeLIB 9.423e-01 4.738e-01 1.989 0.046817 *
## Year_fct2017:StationCodeLIB -3.952e-01 4.620e-01 -0.855 0.392434
## Year_fct2018:StationCodeLIB -2.200e+00 5.051e-01 -4.356 1.36e-05 ***
## Year_fct2019:StationCodeLIB -9.373e-01 5.127e-01 -1.828 0.067617 .
## Year_fct2016:StationCodeRVB 3.293e-01 5.016e-01 0.657 0.511538
## Year_fct2017:StationCodeRVB -5.462e-01 4.829e-01 -1.131 0.258170
## Year_fct2018:StationCodeRVB -4.555e-01 4.624e-01 -0.985 0.324691
## Year_fct2019:StationCodeRVB -6.753e-01 4.697e-01 -1.438 0.150619
## Flow:StationCodeRD22 -4.185e-04 2.542e-04 -1.646 0.099770 .
## Flow:StationCodeI80 -1.249e-03 2.832e-04 -4.412 1.06e-05 ***
## Flow:StationCodeLIS -2.168e-05 2.860e-04 -0.076 0.939582
## Flow:StationCodeSTTD 1.431e-03 2.903e-04 4.930 8.61e-07 ***
## Flow:StationCodeLIB 3.575e-04 1.821e-04 1.964 0.049625 *
## Flow:StationCodeRVB 1.822e-04 1.771e-04 1.029 0.303658
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.675 3.675 11.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.819
## Scale est. = 0.22342 n = 3478
appraise(m_gamm_q_ar2_gam)
k.check(m_gamm_q_ar2_gam)
## k' edf k-index p-value
## s(DOY) 4 3.674585 0.803638 0
draw(m_gamm_q_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_q_ar2)
Box.test(resid_q_ar2, lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: resid_q_ar2
## X-squared = 29.682, df = 20, p-value = 0.07517
The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ (Year_fct + Flow + StationCode)^2 + s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 0.407 0.803
## Flow 1 1.056 0.304
## StationCode 6 14.929 < 2e-16
## Year_fct:Flow 4 0.408 0.803
## Year_fct:StationCode 24 3.933 3.9e-10
## Flow:StationCode 6 16.764 < 2e-16
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.675 3.675 11.3 <2e-16
The Year:StationCode and Flow:StationCode interactions were significant among the interaction terms of the model. Let’s remove the Year:Flow interaction term and check the model summary and diagnostics.
m_gamm_q_ar2b <- gamm(
Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY, k = 5),
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
# Pull out the residuals and the GAM model
resid_q_ar2b <- residuals(m_gamm_q_ar2b$lme, type = "normalized")
m_gamm_q_ar2b_gam <- m_gamm_q_ar2b$gam
summary(m_gamm_q_ar2b_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY,
## k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.180e+00 2.920e-01 31.445 < 2e-16 ***
## Year_fct2016 -3.484e-02 3.758e-01 -0.093 0.926138
## Year_fct2017 3.299e-01 3.584e-01 0.920 0.357500
## Year_fct2018 1.278e-01 3.493e-01 0.366 0.714492
## Year_fct2019 8.237e-02 3.503e-01 0.235 0.814093
## StationCodeRD22 3.144e-01 3.568e-01 0.881 0.378251
## StationCodeI80 7.079e-01 3.669e-01 1.930 0.053736 .
## StationCodeLIS -5.870e-01 3.494e-01 -1.680 0.093030 .
## StationCodeSTTD -4.804e-01 3.574e-01 -1.344 0.178970
## StationCodeLIB -1.239e+00 3.563e-01 -3.477 0.000514 ***
## StationCodeRVB -1.490e+00 3.492e-01 -4.266 2.04e-05 ***
## Flow -1.903e-04 1.769e-04 -1.076 0.282092
## Year_fct2016:StationCodeRD22 -4.388e-01 4.821e-01 -0.910 0.362757
## Year_fct2017:StationCodeRD22 -5.866e-01 4.589e-01 -1.278 0.201177
## Year_fct2018:StationCodeRD22 -4.239e-01 4.483e-01 -0.946 0.344444
## Year_fct2019:StationCodeRD22 -1.649e-01 4.501e-01 -0.366 0.714156
## Year_fct2016:StationCodeI80 1.310e-01 4.893e-01 0.268 0.788936
## Year_fct2017:StationCodeI80 -3.767e-01 4.666e-01 -0.807 0.419558
## Year_fct2018:StationCodeI80 -6.194e-01 4.562e-01 -1.358 0.174616
## Year_fct2019:StationCodeI80 3.978e-01 4.579e-01 0.869 0.385022
## Year_fct2016:StationCodeLIS 6.924e-01 4.697e-01 1.474 0.140536
## Year_fct2017:StationCodeLIS 1.681e-01 4.563e-01 0.368 0.712634
## Year_fct2018:StationCodeLIS -1.545e-01 4.442e-01 -0.348 0.728056
## Year_fct2019:StationCodeLIS 4.184e-01 4.493e-01 0.931 0.351866
## Year_fct2016:StationCodeSTTD 3.480e-01 4.846e-01 0.718 0.472790
## Year_fct2017:StationCodeSTTD -5.708e-01 4.911e-01 -1.162 0.245159
## Year_fct2018:StationCodeSTTD -8.680e-01 4.657e-01 -1.864 0.062400 .
## Year_fct2019:StationCodeSTTD -1.191e+00 4.581e-01 -2.600 0.009364 **
## Year_fct2016:StationCodeLIB 9.536e-01 4.727e-01 2.017 0.043741 *
## Year_fct2017:StationCodeLIB -3.846e-01 4.611e-01 -0.834 0.404321
## Year_fct2018:StationCodeLIB -2.191e+00 5.042e-01 -4.345 1.43e-05 ***
## Year_fct2019:StationCodeLIB -9.096e-01 5.118e-01 -1.777 0.075641 .
## Year_fct2016:StationCodeRVB 3.114e-01 4.819e-01 0.646 0.518156
## Year_fct2017:StationCodeRVB -5.575e-01 4.532e-01 -1.230 0.218786
## Year_fct2018:StationCodeRVB -4.597e-01 4.419e-01 -1.040 0.298248
## Year_fct2019:StationCodeRVB -8.461e-01 4.498e-01 -1.881 0.060059 .
## StationCodeRD22:Flow -4.181e-04 2.541e-04 -1.645 0.100007
## StationCodeI80:Flow -1.252e-03 2.831e-04 -4.424 1.00e-05 ***
## StationCodeLIS:Flow -2.371e-05 2.859e-04 -0.083 0.933893
## StationCodeSTTD:Flow 1.429e-03 2.902e-04 4.924 8.88e-07 ***
## StationCodeLIB:Flow 3.608e-04 1.819e-04 1.983 0.047421 *
## StationCodeRVB:Flow 1.840e-04 1.770e-04 1.039 0.298683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.68 3.68 11.47 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.819
## Scale est. = 0.22326 n = 3478
appraise(m_gamm_q_ar2b_gam)
k.check(m_gamm_q_ar2b_gam)
## k' edf k-index p-value
## s(DOY) 4 3.679956 0.8026532 0
draw(m_gamm_q_ar2b_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_ar2b_gam, pages = 1, all.terms = TRUE)
acf(resid_q_ar2b)
Box.test(resid_q_ar2b, lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: resid_q_ar2b
## X-squared = 29.646, df = 20, p-value = 0.0758
anova(m_gamm_q_ar2b_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + s(DOY,
## k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 0.406 0.805
## StationCode 6 14.873 < 2e-16
## Flow 1 1.157 0.282
## Year_fct:StationCode 24 4.102 8.4e-11
## StationCode:Flow 6 16.865 < 2e-16
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.68 3.68 11.47 <2e-16
This model still looks pretty good. We’ll use this one for the model selection process.
Before we add both flow and water temperature as continuous predictors in the model, let’s check if they are correlated with each other.
df_chla_wt_q %>%
ggplot(aes(x = Flow, y = WaterTemp)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(vars(StationCode), scales = "free") +
theme_bw()
Flow and water temperature appear to have no correlation.
We’ll try running a GAM using the model structure from the flow model and adding two-way interactions between Year:WaterTemp and Station:WaterTemp. First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_q_wt <- gam(
Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)
Lets look at the model summary and diagnostics:
summary(m_gam_q_wt)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct *
## WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.2662032 0.2178374 47.128 < 2e-16 ***
## Year_fct2016 -0.6484169 0.3059264 -2.120 0.034118 *
## Year_fct2017 -0.0609459 0.1748727 -0.349 0.727474
## Year_fct2018 -0.1613266 0.1680842 -0.960 0.337226
## Year_fct2019 0.5338094 0.1622929 3.289 0.001015 **
## StationCodeRD22 -1.0478545 0.1618085 -6.476 1.08e-10 ***
## StationCodeI80 -0.3198595 0.1552094 -2.061 0.039395 *
## StationCodeLIS -2.5319716 0.1616506 -15.663 < 2e-16 ***
## StationCodeSTTD -1.7388040 0.1881330 -9.242 < 2e-16 ***
## StationCodeLIB -3.9184097 0.2367461 -16.551 < 2e-16 ***
## StationCodeRVB -3.4749646 0.1941132 -17.902 < 2e-16 ***
## Flow -0.0011256 0.0001014 -11.096 < 2e-16 ***
## WaterTemp -0.0433571 0.0092626 -4.681 2.97e-06 ***
## Year_fct2016:StationCodeRD22 -0.5347596 0.0910571 -5.873 4.69e-09 ***
## Year_fct2017:StationCodeRD22 -0.6200951 0.0889603 -6.970 3.77e-12 ***
## Year_fct2018:StationCodeRD22 -0.3517822 0.0834277 -4.217 2.54e-05 ***
## Year_fct2019:StationCodeRD22 -0.1482813 0.0833741 -1.779 0.075410 .
## Year_fct2016:StationCodeI80 0.2012539 0.0934108 2.155 0.031270 *
## Year_fct2017:StationCodeI80 -0.2400564 0.0912588 -2.631 0.008564 **
## Year_fct2018:StationCodeI80 -0.4104957 0.0855122 -4.800 1.65e-06 ***
## Year_fct2019:StationCodeI80 0.5658673 0.0856073 6.610 4.44e-11 ***
## Year_fct2016:StationCodeLIS 0.3694915 0.0892519 4.140 3.56e-05 ***
## Year_fct2017:StationCodeLIS 0.2490534 0.0878966 2.833 0.004631 **
## Year_fct2018:StationCodeLIS -0.0695068 0.0825594 -0.842 0.399902
## Year_fct2019:StationCodeLIS 0.5781116 0.0838872 6.892 6.54e-12 ***
## Year_fct2016:StationCodeSTTD 0.2714631 0.0931049 2.916 0.003572 **
## Year_fct2017:StationCodeSTTD -0.1673873 0.0975732 -1.716 0.086343 .
## Year_fct2018:StationCodeSTTD -0.6929241 0.0871555 -7.950 2.50e-15 ***
## Year_fct2019:StationCodeSTTD -1.1345661 0.0858236 -13.220 < 2e-16 ***
## Year_fct2016:StationCodeLIB 0.8903568 0.0949408 9.378 < 2e-16 ***
## Year_fct2017:StationCodeLIB -0.2315049 0.0969164 -2.389 0.016962 *
## Year_fct2018:StationCodeLIB -2.0080384 0.1150845 -17.448 < 2e-16 ***
## Year_fct2019:StationCodeLIB -0.8582224 0.1196082 -7.175 8.81e-13 ***
## Year_fct2016:StationCodeRVB 0.1597391 0.1027829 1.554 0.120243
## Year_fct2017:StationCodeRVB -0.3753724 0.1038743 -3.614 0.000306 ***
## Year_fct2018:StationCodeRVB -0.2920268 0.0889282 -3.284 0.001034 **
## Year_fct2019:StationCodeRVB -0.6540161 0.0928199 -7.046 2.22e-12 ***
## StationCodeRD22:Flow -0.0003430 0.0001340 -2.560 0.010519 *
## StationCodeI80:Flow -0.0008049 0.0001352 -5.951 2.93e-09 ***
## StationCodeLIS:Flow 0.0012604 0.0001350 9.336 < 2e-16 ***
## StationCodeSTTD:Flow 0.0030649 0.0001386 22.121 < 2e-16 ***
## StationCodeLIB:Flow 0.0013482 0.0001236 10.908 < 2e-16 ***
## StationCodeRVB:Flow 0.0011204 0.0001014 11.046 < 2e-16 ***
## Year_fct2016:WaterTemp 0.0278927 0.0123785 2.253 0.024302 *
## Year_fct2017:WaterTemp 0.0125116 0.0069925 1.789 0.073658 .
## Year_fct2018:WaterTemp 0.0103364 0.0068995 1.498 0.134191
## Year_fct2019:WaterTemp -0.0253859 0.0065403 -3.881 0.000106 ***
## StationCodeRD22:WaterTemp 0.0637778 0.0067035 9.514 < 2e-16 ***
## StationCodeI80:WaterTemp 0.0389816 0.0064179 6.074 1.39e-09 ***
## StationCodeLIS:WaterTemp 0.0857066 0.0066872 12.816 < 2e-16 ***
## StationCodeSTTD:WaterTemp 0.0496447 0.0081526 6.089 1.26e-09 ***
## StationCodeLIB:WaterTemp 0.1230615 0.0120862 10.182 < 2e-16 ***
## StationCodeRVB:WaterTemp 0.0833896 0.0086653 9.623 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.781 3.974 53.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.864 Deviance explained = 86.6%
## -REML = 1615.8 Scale est. = 0.13435 n = 3478
appraise(m_gam_q_wt)
k.check(m_gam_q_wt)
## k' edf k-index p-value
## s(DOY) 4 3.78136 0.9264289 0
draw(m_gam_q_wt, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_q_wt, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_q_wt))
Box.test(residuals(m_gam_q_wt), lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: residuals(m_gam_q_wt)
## X-squared = 8803.9, df = 20, p-value < 2.2e-16
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_q_wt <- as.formula(
"Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)"
)
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_wt <- gamm(
f_gam_q_wt,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_q_wt_ar1 <- gamm(
f_gam_q_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_wt_ar2 <- gamm(
f_gam_q_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_wt_ar3 <- gamm(
f_gam_q_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_q_wt$lme,
m_gamm_q_wt_ar1$lme,
m_gamm_q_wt_ar2$lme,
m_gamm_q_wt_ar3$lme
)
## Model df AIC BIC logLik Test L.Ratio
## m_gamm_q_wt$lme 1 56 3343.533 3687.292 -1615.766
## m_gamm_q_wt_ar1$lme 2 57 -2152.733 -1802.835 1133.367 1 vs 2 5498.266
## m_gamm_q_wt_ar2$lme 3 58 -2174.391 -1818.354 1145.195 2 vs 3 23.657
## m_gamm_q_wt_ar3$lme 4 59 -2173.370 -1811.195 1145.685 3 vs 4 0.980
## p-value
## m_gamm_q_wt$lme
## m_gamm_q_wt_ar1$lme <.0001
## m_gamm_q_wt_ar2$lme <.0001
## m_gamm_q_wt_ar3$lme 0.3222
It looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
# Pull out the residuals and the GAM model
resid_q_wt_ar2 <- residuals(m_gamm_q_wt_ar2$lme, type = "normalized")
m_gamm_q_wt_ar2_gam <- m_gamm_q_wt_ar2$gam
summary(m_gamm_q_wt_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct *
## WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.894e+00 5.959e-01 14.925 < 2e-16 ***
## Year_fct2016 1.423e+00 6.306e-01 2.256 0.024113 *
## Year_fct2017 4.921e-02 5.714e-01 0.086 0.931374
## Year_fct2018 -3.647e-01 5.778e-01 -0.631 0.527977
## Year_fct2019 -2.601e-02 5.490e-01 -0.047 0.962218
## StationCodeRD22 -1.128e+00 6.141e-01 -1.838 0.066189 .
## StationCodeI80 8.206e-01 5.957e-01 1.378 0.168447
## StationCodeLIS -1.605e+00 6.037e-01 -2.659 0.007879 **
## StationCodeSTTD -1.252e+00 6.252e-01 -2.003 0.045228 *
## StationCodeLIB -9.809e-01 6.957e-01 -1.410 0.158668
## StationCodeRVB -2.142e+00 7.531e-01 -2.845 0.004470 **
## Flow -1.699e-04 1.767e-04 -0.961 0.336414
## WaterTemp 1.399e-02 2.161e-02 0.648 0.517248
## Year_fct2016:StationCodeRD22 -6.003e-01 5.380e-01 -1.116 0.264600
## Year_fct2017:StationCodeRD22 -6.123e-01 5.116e-01 -1.197 0.231484
## Year_fct2018:StationCodeRD22 -3.207e-01 5.017e-01 -0.639 0.522705
## Year_fct2019:StationCodeRD22 -1.397e-01 5.030e-01 -0.278 0.781210
## Year_fct2016:StationCodeI80 2.047e-01 5.478e-01 0.374 0.708640
## Year_fct2017:StationCodeI80 -2.886e-01 5.216e-01 -0.553 0.580145
## Year_fct2018:StationCodeI80 -5.347e-01 5.115e-01 -1.045 0.295920
## Year_fct2019:StationCodeI80 4.585e-01 5.129e-01 0.894 0.371361
## Year_fct2016:StationCodeLIS 5.821e-01 5.244e-01 1.110 0.267053
## Year_fct2017:StationCodeLIS 1.927e-01 5.089e-01 0.379 0.704890
## Year_fct2018:StationCodeLIS -1.008e-01 4.971e-01 -0.203 0.839336
## Year_fct2019:StationCodeLIS 4.263e-01 5.023e-01 0.849 0.396085
## Year_fct2016:StationCodeSTTD 2.585e-01 5.418e-01 0.477 0.633378
## Year_fct2017:StationCodeSTTD -5.788e-01 5.483e-01 -1.056 0.291233
## Year_fct2018:StationCodeSTTD -8.371e-01 5.208e-01 -1.607 0.108099
## Year_fct2019:StationCodeSTTD -1.172e+00 5.124e-01 -2.288 0.022224 *
## Year_fct2016:StationCodeLIB 8.689e-01 5.287e-01 1.644 0.100355
## Year_fct2017:StationCodeLIB -2.787e-01 5.147e-01 -0.541 0.588263
## Year_fct2018:StationCodeLIB -2.096e+00 5.618e-01 -3.731 0.000194 ***
## Year_fct2019:StationCodeLIB -8.120e-01 5.685e-01 -1.428 0.153286
## Year_fct2016:StationCodeRVB 1.962e-01 5.384e-01 0.364 0.715519
## Year_fct2017:StationCodeRVB -4.613e-01 5.059e-01 -0.912 0.361903
## Year_fct2018:StationCodeRVB -3.546e-01 4.956e-01 -0.715 0.474415
## Year_fct2019:StationCodeRVB -7.574e-01 5.041e-01 -1.503 0.133011
## StationCodeRD22:Flow -4.412e-04 2.550e-04 -1.730 0.083665 .
## StationCodeI80:Flow -1.205e-03 2.852e-04 -4.226 2.44e-05 ***
## StationCodeLIS:Flow -8.722e-06 2.877e-04 -0.030 0.975817
## StationCodeSTTD:Flow 1.426e-03 2.923e-04 4.878 1.12e-06 ***
## StationCodeLIB:Flow 3.233e-04 1.821e-04 1.775 0.075912 .
## StationCodeRVB:Flow 1.655e-04 1.768e-04 0.936 0.349418
## Year_fct2016:WaterTemp -5.896e-02 2.002e-02 -2.945 0.003250 **
## Year_fct2017:WaterTemp 1.079e-02 1.791e-02 0.603 0.546631
## Year_fct2018:WaterTemp 2.128e-02 1.913e-02 1.113 0.265978
## Year_fct2019:WaterTemp 3.786e-03 1.708e-02 0.222 0.824591
## StationCodeRD22:WaterTemp 6.651e-02 2.099e-02 3.169 0.001541 **
## StationCodeI80:WaterTemp -8.905e-03 1.954e-02 -0.456 0.648553
## StationCodeLIS:WaterTemp 4.642e-02 2.053e-02 2.261 0.023804 *
## StationCodeSTTD:WaterTemp 3.510e-02 2.186e-02 1.606 0.108445
## StationCodeLIB:WaterTemp -1.635e-02 2.705e-02 -0.605 0.545485
## StationCodeRVB:WaterTemp 2.870e-02 2.960e-02 0.970 0.332241
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.664 3.664 8.628 2.2e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.807
## Scale est. = 0.2515 n = 3478
appraise(m_gamm_q_wt_ar2_gam)
k.check(m_gamm_q_wt_ar2_gam)
## k' edf k-index p-value
## s(DOY) 4 3.663875 0.7855439 0
draw(m_gamm_q_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_wt_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_q_wt_ar2)
Box.test(resid_q_wt_ar2, lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: resid_q_wt_ar2
## X-squared = 30.519, df = 20, p-value = 0.06187
The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_wt_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode + Year_fct *
## WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 2.531 0.038608
## StationCode 6 6.652 5.20e-07
## Flow 1 0.924 0.336414
## WaterTemp 1 0.419 0.517248
## Year_fct:StationCode 24 3.165 3.27e-07
## StationCode:Flow 6 15.268 < 2e-16
## Year_fct:WaterTemp 4 4.659 0.000945
## StationCode:WaterTemp 6 4.669 9.75e-05
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.664 3.664 8.628 2.2e-06
All interaction terms are significant in this model, so we’ll use this one for the model selection process.
Let’s look at the relationship between chlorophyll and flow for each station-flow action period combination.
df_chla_wt_q %>%
ggplot(aes(x = Flow, y = Chla)) +
geom_point() +
geom_smooth(formula = "y ~ x") +
facet_wrap(
vars(StationCode, FlowActionPeriod),
ncol = 3,
scales = "free",
labeller = labeller(.multi_line = FALSE)
) +
theme_bw()
It’s possible that the relationship between chlorophyll and Flow differs between Flow Action Periods at each station. We’ll try to account for this by running a GAM including a two-way interaction between Year and Station, a three-way interaction between Flow, Station, and Flow Action Period, and a smooth term for day of year to account for seasonality (restricting the k-value to 5 to reduce overfitting). We won’t include the interaction between Year and Flow Action Period since it was barely significant in the categorical model and we’re primarily including the Flow Action Period to allow for different slopes for flow at each station. First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_q_fa <- gam(
Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)
Lets look at the model summary and diagnostics:
summary(m_gam_q_fa)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 9.0200439 0.0839467 107.450
## Year_fct2016 -0.0752673 0.0817958 -0.920
## Year_fct2017 0.1772134 0.0707049 2.506
## Year_fct2018 -0.0051101 0.0673845 -0.076
## Year_fct2019 -0.0504808 0.0701749 -0.719
## StationCodeRD22 0.6048920 0.1021743 5.920
## StationCodeI80 0.7075466 0.1053335 6.717
## StationCodeLIS -0.2534964 0.1008109 -2.515
## StationCodeSTTD -0.4888704 0.1077499 -4.537
## StationCodeLIB -0.5233575 0.1654118 -3.164
## StationCodeRVB -1.3296625 0.1140411 -11.660
## Flow 0.0003845 0.0006336 0.607
## FlowActionPeriodDuring 0.2206777 0.0876376 2.518
## FlowActionPeriodAfter 0.8225799 0.0753879 10.911
## Year_fct2016:StationCodeRD22 -0.3499073 0.0955876 -3.661
## Year_fct2017:StationCodeRD22 -0.5693333 0.0856991 -6.643
## Year_fct2018:StationCodeRD22 -0.3195979 0.0814168 -3.925
## Year_fct2019:StationCodeRD22 -0.0887040 0.0850886 -1.042
## Year_fct2016:StationCodeI80 0.3049614 0.0984091 3.099
## Year_fct2017:StationCodeI80 -0.2386526 0.0866819 -2.753
## Year_fct2018:StationCodeI80 -0.3895701 0.0840942 -4.633
## Year_fct2019:StationCodeI80 0.6211290 0.0862887 7.198
## Year_fct2016:StationCodeLIS 0.3767691 0.0932856 4.039
## Year_fct2017:StationCodeLIS 0.2509619 0.0859328 2.920
## Year_fct2018:StationCodeLIS 0.0062195 0.0810131 0.077
## Year_fct2019:StationCodeLIS 0.6367041 0.0862458 7.382
## Year_fct2016:StationCodeSTTD 0.2375629 0.0960514 2.473
## Year_fct2017:StationCodeSTTD -0.0453964 0.0976543 -0.465
## Year_fct2018:StationCodeSTTD -0.5512862 0.0857425 -6.430
## Year_fct2019:StationCodeSTTD -1.0150911 0.0885178 -11.468
## Year_fct2016:StationCodeLIB 0.8861875 0.0923868 9.592
## Year_fct2017:StationCodeLIB -0.2709091 0.0930803 -2.910
## Year_fct2018:StationCodeLIB -2.0402939 0.1103463 -18.490
## Year_fct2019:StationCodeLIB -1.0231592 0.1184799 -8.636
## Year_fct2016:StationCodeRVB 0.1576659 0.1008718 1.563
## Year_fct2017:StationCodeRVB -0.4120981 0.0981257 -4.200
## Year_fct2018:StationCodeRVB -0.2930854 0.0862916 -3.396
## Year_fct2019:StationCodeRVB -0.6454502 0.0914059 -7.061
## StationCodeRD22:Flow 0.0003861 0.0014058 0.275
## StationCodeI80:Flow 0.0018616 0.0018957 0.982
## StationCodeLIS:Flow 0.0030944 0.0009665 3.202
## StationCodeSTTD:Flow 0.0023417 0.0011000 2.129
## StationCodeLIB:Flow 0.0001091 0.0006393 0.171
## StationCodeRVB:Flow -0.0003875 0.0006337 -0.612
## Flow:FlowActionPeriodDuring -0.0013504 0.0006478 -2.084
## Flow:FlowActionPeriodAfter -0.0083453 0.0011337 -7.361
## StationCodeRD22:FlowActionPeriodDuring -0.2621276 0.1125150 -2.330
## StationCodeI80:FlowActionPeriodDuring -0.1046615 0.1170017 -0.895
## StationCodeLIS:FlowActionPeriodDuring 0.3486851 0.1118607 3.117
## StationCodeSTTD:FlowActionPeriodDuring 0.4101368 0.1221588 3.357
## StationCodeLIB:FlowActionPeriodDuring -0.3116617 0.1631523 -1.910
## StationCodeRVB:FlowActionPeriodDuring -0.2557270 0.1403998 -1.821
## StationCodeRD22:FlowActionPeriodAfter -0.6738249 0.0876828 -7.685
## StationCodeI80:FlowActionPeriodAfter -0.4696770 0.0926079 -5.072
## StationCodeLIS:FlowActionPeriodAfter -0.6118818 0.0862036 -7.098
## StationCodeSTTD:FlowActionPeriodAfter -0.4414052 0.0972779 -4.538
## StationCodeLIB:FlowActionPeriodAfter -1.9012296 0.1443312 -13.173
## StationCodeRVB:FlowActionPeriodAfter -0.9357160 0.1208432 -7.743
## StationCodeRD22:Flow:FlowActionPeriodDuring -0.0007004 0.0014165 -0.494
## StationCodeI80:Flow:FlowActionPeriodDuring -0.0028106 0.0019000 -1.479
## StationCodeLIS:Flow:FlowActionPeriodDuring -0.0035976 0.0009869 -3.645
## StationCodeSTTD:Flow:FlowActionPeriodDuring -0.0007149 0.0011208 -0.638
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.0011745 0.0006559 1.791
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.0013536 0.0006479 2.089
## StationCodeRD22:Flow:FlowActionPeriodAfter -0.0017581 0.0023433 -0.750
## StationCodeI80:Flow:FlowActionPeriodAfter -0.0055463 0.0025798 -2.150
## StationCodeLIS:Flow:FlowActionPeriodAfter 0.0125081 0.0015103 8.282
## StationCodeSTTD:Flow:FlowActionPeriodAfter 0.0118661 0.0016282 7.288
## StationCodeLIB:Flow:FlowActionPeriodAfter 0.0076975 0.0011391 6.758
## StationCodeRVB:Flow:FlowActionPeriodAfter 0.0083444 0.0011336 7.361
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2016 0.357541
## Year_fct2017 0.012244 *
## Year_fct2018 0.939555
## Year_fct2019 0.471970
## StationCodeRD22 3.53e-09 ***
## StationCodeI80 2.16e-11 ***
## StationCodeLIS 0.011963 *
## StationCodeSTTD 5.90e-06 ***
## StationCodeLIB 0.001570 **
## StationCodeRVB < 2e-16 ***
## Flow 0.543936
## FlowActionPeriodDuring 0.011846 *
## FlowActionPeriodAfter < 2e-16 ***
## Year_fct2016:StationCodeRD22 0.000255 ***
## Year_fct2017:StationCodeRD22 3.55e-11 ***
## Year_fct2018:StationCodeRD22 8.83e-05 ***
## Year_fct2019:StationCodeRD22 0.297259
## Year_fct2016:StationCodeI80 0.001958 **
## Year_fct2017:StationCodeI80 0.005933 **
## Year_fct2018:StationCodeI80 3.75e-06 ***
## Year_fct2019:StationCodeI80 7.47e-13 ***
## Year_fct2016:StationCodeLIS 5.49e-05 ***
## Year_fct2017:StationCodeLIS 0.003518 **
## Year_fct2018:StationCodeLIS 0.938810
## Year_fct2019:StationCodeLIS 1.94e-13 ***
## Year_fct2016:StationCodeSTTD 0.013436 *
## Year_fct2017:StationCodeSTTD 0.642055
## Year_fct2018:StationCodeSTTD 1.46e-10 ***
## Year_fct2019:StationCodeSTTD < 2e-16 ***
## Year_fct2016:StationCodeLIB < 2e-16 ***
## Year_fct2017:StationCodeLIB 0.003632 **
## Year_fct2018:StationCodeLIB < 2e-16 ***
## Year_fct2019:StationCodeLIB < 2e-16 ***
## Year_fct2016:StationCodeRVB 0.118138
## Year_fct2017:StationCodeRVB 2.74e-05 ***
## Year_fct2018:StationCodeRVB 0.000690 ***
## Year_fct2019:StationCodeRVB 1.99e-12 ***
## StationCodeRD22:Flow 0.783602
## StationCodeI80:Flow 0.326154
## StationCodeLIS:Flow 0.001378 **
## StationCodeSTTD:Flow 0.033342 *
## StationCodeLIB:Flow 0.864475
## StationCodeRVB:Flow 0.540909
## Flow:FlowActionPeriodDuring 0.037198 *
## Flow:FlowActionPeriodAfter 2.27e-13 ***
## StationCodeRD22:FlowActionPeriodDuring 0.019880 *
## StationCodeI80:FlowActionPeriodDuring 0.371102
## StationCodeLIS:FlowActionPeriodDuring 0.001841 **
## StationCodeSTTD:FlowActionPeriodDuring 0.000795 ***
## StationCodeLIB:FlowActionPeriodDuring 0.056185 .
## StationCodeRVB:FlowActionPeriodDuring 0.068631 .
## StationCodeRD22:FlowActionPeriodAfter 1.99e-14 ***
## StationCodeI80:FlowActionPeriodAfter 4.15e-07 ***
## StationCodeLIS:FlowActionPeriodAfter 1.53e-12 ***
## StationCodeSTTD:FlowActionPeriodAfter 5.89e-06 ***
## StationCodeLIB:FlowActionPeriodAfter < 2e-16 ***
## StationCodeRVB:FlowActionPeriodAfter 1.27e-14 ***
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.621005
## StationCodeI80:Flow:FlowActionPeriodDuring 0.139151
## StationCodeLIS:Flow:FlowActionPeriodDuring 0.000271 ***
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.523617
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.073440 .
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.036769 *
## StationCodeRD22:Flow:FlowActionPeriodAfter 0.453163
## StationCodeI80:Flow:FlowActionPeriodAfter 0.031634 *
## StationCodeLIS:Flow:FlowActionPeriodAfter < 2e-16 ***
## StationCodeSTTD:Flow:FlowActionPeriodAfter 3.90e-13 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter 1.64e-11 ***
## StationCodeRVB:Flow:FlowActionPeriodAfter 2.28e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.863 3.989 47.56 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.88 Deviance explained = 88.2%
## -REML = 1476.1 Scale est. = 0.11868 n = 3478
appraise(m_gam_q_fa)
k.check(m_gam_q_fa)
## k' edf k-index p-value
## s(DOY) 4 3.862781 0.9687532 0.0375
draw(m_gam_q_fa, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_q_fa, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_q_fa))
Box.test(residuals(m_gam_q_fa), lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: residuals(m_gam_q_fa)
## X-squared = 7116.1, df = 20, p-value < 2.2e-16
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_q_fa <- as.formula("Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + s(DOY, k = 5)")
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_fa <- gamm(
f_gam_q_fa,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_q_fa_ar1 <- gamm(
f_gam_q_fa,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_fa_ar2 <- gamm(
f_gam_q_fa,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_fa_ar3 <- gamm(
f_gam_q_fa,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_q_fa$lme,
m_gamm_q_fa_ar1$lme,
m_gamm_q_fa_ar2$lme,
m_gamm_q_fa_ar3$lme
)
## Model df AIC BIC logLik Test L.Ratio
## m_gamm_q_fa$lme 1 73 3098.134 3545.886 -1476.067
## m_gamm_q_fa_ar1$lme 2 74 -1985.717 -1531.832 1066.859 1 vs 2 5085.851
## m_gamm_q_fa_ar2$lme 3 75 -2005.949 -1545.930 1077.975 2 vs 3 22.232
## m_gamm_q_fa_ar3$lme 4 76 -2004.214 -1538.061 1078.107 3 vs 4 0.265
## p-value
## m_gamm_q_fa$lme
## m_gamm_q_fa_ar1$lme <.0001
## m_gamm_q_fa_ar2$lme <.0001
## m_gamm_q_fa_ar3$lme 0.6067
Again, it looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
# Pull out the residuals and the GAM model
resid_q_fa_ar2 <- residuals(m_gamm_q_fa_ar2$lme, type = "normalized")
m_gamm_q_fa_ar2_gam <- m_gamm_q_fa_ar2$gam
summary(m_gamm_q_fa_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 9.140e+00 3.300e-01 27.694
## Year_fct2016 -8.591e-02 4.210e-01 -0.204
## Year_fct2017 3.071e-01 3.994e-01 0.769
## Year_fct2018 1.029e-01 3.897e-01 0.264
## Year_fct2019 6.427e-02 3.912e-01 0.164
## StationCodeRD22 2.084e-01 4.107e-01 0.508
## StationCodeI80 3.448e-01 4.231e-01 0.815
## StationCodeLIS -6.571e-01 3.993e-01 -1.646
## StationCodeSTTD -2.267e-01 4.099e-01 -0.553
## StationCodeLIB -9.494e-01 4.148e-01 -2.289
## StationCodeRVB -1.386e+00 4.029e-01 -3.441
## Flow -1.686e-04 3.656e-04 -0.461
## FlowActionPeriodDuring 8.981e-02 1.096e-01 0.819
## FlowActionPeriodAfter 1.447e-01 1.313e-01 1.102
## Year_fct2016:StationCodeRD22 -4.345e-01 5.402e-01 -0.804
## Year_fct2017:StationCodeRD22 -5.513e-01 5.127e-01 -1.075
## Year_fct2018:StationCodeRD22 -3.939e-01 5.015e-01 -0.785
## Year_fct2019:StationCodeRD22 -1.286e-01 5.039e-01 -0.255
## Year_fct2016:StationCodeI80 1.236e-01 5.486e-01 0.225
## Year_fct2017:StationCodeI80 -3.138e-01 5.212e-01 -0.602
## Year_fct2018:StationCodeI80 -5.976e-01 5.102e-01 -1.171
## Year_fct2019:StationCodeI80 4.697e-01 5.125e-01 0.916
## Year_fct2016:StationCodeLIS 7.761e-01 5.260e-01 1.475
## Year_fct2017:StationCodeLIS 2.046e-01 5.104e-01 0.401
## Year_fct2018:StationCodeLIS -1.434e-01 4.971e-01 -0.288
## Year_fct2019:StationCodeLIS 4.044e-01 5.029e-01 0.804
## Year_fct2016:StationCodeSTTD 3.971e-01 5.426e-01 0.732
## Year_fct2017:StationCodeSTTD -6.679e-01 5.490e-01 -1.217
## Year_fct2018:StationCodeSTTD -9.501e-01 5.208e-01 -1.824
## Year_fct2019:StationCodeSTTD -1.144e+00 5.126e-01 -2.232
## Year_fct2016:StationCodeLIB 9.970e-01 5.295e-01 1.883
## Year_fct2017:StationCodeLIB -4.086e-01 5.154e-01 -0.793
## Year_fct2018:StationCodeLIB -2.176e+00 5.607e-01 -3.881
## Year_fct2019:StationCodeLIB -9.519e-01 5.694e-01 -1.672
## Year_fct2016:StationCodeRVB 3.678e-01 5.390e-01 0.682
## Year_fct2017:StationCodeRVB -5.412e-01 5.064e-01 -1.069
## Year_fct2018:StationCodeRVB -4.418e-01 4.948e-01 -0.893
## Year_fct2019:StationCodeRVB -8.452e-01 5.035e-01 -1.678
## StationCodeRD22:Flow 6.557e-04 1.189e-03 0.552
## StationCodeI80:Flow 4.714e-03 1.602e-03 2.943
## StationCodeLIS:Flow -8.155e-04 8.945e-04 -0.912
## StationCodeSTTD:Flow -1.186e-04 1.026e-03 -0.116
## StationCodeLIB:Flow 4.674e-04 3.725e-04 1.255
## StationCodeRVB:Flow 1.588e-04 3.657e-04 0.434
## Flow:FlowActionPeriodDuring -5.381e-05 3.758e-04 -0.143
## Flow:FlowActionPeriodAfter 7.104e-05 7.263e-04 0.098
## StationCodeRD22:FlowActionPeriodDuring 3.777e-02 1.509e-01 0.250
## StationCodeI80:FlowActionPeriodDuring 3.444e-01 1.540e-01 2.236
## StationCodeLIS:FlowActionPeriodDuring -1.107e-01 1.328e-01 -0.834
## StationCodeSTTD:FlowActionPeriodDuring 7.934e-03 1.340e-01 0.059
## StationCodeLIB:FlowActionPeriodDuring -4.421e-01 1.783e-01 -2.480
## StationCodeRVB:FlowActionPeriodDuring -1.782e-01 1.846e-01 -0.965
## StationCodeRD22:FlowActionPeriodAfter 8.258e-02 1.954e-01 0.423
## StationCodeI80:FlowActionPeriodAfter 4.533e-01 1.984e-01 2.285
## StationCodeLIS:FlowActionPeriodAfter 1.976e-02 1.635e-01 0.121
## StationCodeSTTD:FlowActionPeriodAfter -4.926e-01 1.658e-01 -2.972
## StationCodeLIB:FlowActionPeriodAfter -4.578e-01 1.974e-01 -2.320
## StationCodeRVB:FlowActionPeriodAfter -2.304e-01 1.999e-01 -1.153
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.040e-03 1.172e-03 -0.887
## StationCodeI80:Flow:FlowActionPeriodDuring -5.914e-03 1.569e-03 -3.771
## StationCodeLIS:Flow:FlowActionPeriodDuring 1.118e-03 9.342e-04 1.197
## StationCodeSTTD:Flow:FlowActionPeriodDuring 9.933e-04 1.061e-03 0.936
## StationCodeLIB:Flow:FlowActionPeriodDuring -1.879e-04 3.889e-04 -0.483
## StationCodeRVB:Flow:FlowActionPeriodDuring 6.217e-05 3.761e-04 0.165
## StationCodeRD22:Flow:FlowActionPeriodAfter -4.517e-05 2.549e-03 -0.018
## StationCodeI80:Flow:FlowActionPeriodAfter -8.338e-03 2.762e-03 -3.019
## StationCodeLIS:Flow:FlowActionPeriodAfter -4.819e-05 1.269e-03 -0.038
## StationCodeSTTD:Flow:FlowActionPeriodAfter 4.786e-03 1.381e-03 3.465
## StationCodeLIB:Flow:FlowActionPeriodAfter -2.209e-04 7.319e-04 -0.302
## StationCodeRVB:Flow:FlowActionPeriodAfter -6.675e-05 7.263e-04 -0.092
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2016 0.838299
## Year_fct2017 0.441995
## Year_fct2018 0.791690
## Year_fct2019 0.869520
## StationCodeRD22 0.611816
## StationCodeI80 0.415252
## StationCodeLIS 0.099916 .
## StationCodeSTTD 0.580352
## StationCodeLIB 0.022136 *
## StationCodeRVB 0.000587 ***
## Flow 0.644739
## FlowActionPeriodDuring 0.412618
## FlowActionPeriodAfter 0.270588
## Year_fct2016:StationCodeRD22 0.421258
## Year_fct2017:StationCodeRD22 0.282374
## Year_fct2018:StationCodeRD22 0.432237
## Year_fct2019:StationCodeRD22 0.798634
## Year_fct2016:StationCodeI80 0.821794
## Year_fct2017:StationCodeI80 0.547132
## Year_fct2018:StationCodeI80 0.241558
## Year_fct2019:StationCodeI80 0.359519
## Year_fct2016:StationCodeLIS 0.140224
## Year_fct2017:StationCodeLIS 0.688545
## Year_fct2018:StationCodeLIS 0.772996
## Year_fct2019:StationCodeLIS 0.421423
## Year_fct2016:StationCodeSTTD 0.464298
## Year_fct2017:StationCodeSTTD 0.223845
## Year_fct2018:StationCodeSTTD 0.068178 .
## Year_fct2019:StationCodeSTTD 0.025692 *
## Year_fct2016:StationCodeLIB 0.059765 .
## Year_fct2017:StationCodeLIB 0.427980
## Year_fct2018:StationCodeLIB 0.000106 ***
## Year_fct2019:StationCodeLIB 0.094640 .
## Year_fct2016:StationCodeRVB 0.495041
## Year_fct2017:StationCodeRVB 0.285267
## Year_fct2018:StationCodeRVB 0.371985
## Year_fct2019:StationCodeRVB 0.093350 .
## StationCodeRD22:Flow 0.581215
## StationCodeI80:Flow 0.003268 **
## StationCodeLIS:Flow 0.362050
## StationCodeSTTD:Flow 0.908046
## StationCodeLIB:Flow 0.209715
## StationCodeRVB:Flow 0.664224
## Flow:FlowActionPeriodDuring 0.886142
## Flow:FlowActionPeriodAfter 0.922079
## StationCodeRD22:FlowActionPeriodDuring 0.802365
## StationCodeI80:FlowActionPeriodDuring 0.025409 *
## StationCodeLIS:FlowActionPeriodDuring 0.404406
## StationCodeSTTD:FlowActionPeriodDuring 0.952780
## StationCodeLIB:FlowActionPeriodDuring 0.013178 *
## StationCodeRVB:FlowActionPeriodDuring 0.334379
## StationCodeRD22:FlowActionPeriodAfter 0.672656
## StationCodeI80:FlowActionPeriodAfter 0.022380 *
## StationCodeLIS:FlowActionPeriodAfter 0.903816
## StationCodeSTTD:FlowActionPeriodAfter 0.002984 **
## StationCodeLIB:FlowActionPeriodAfter 0.020421 *
## StationCodeRVB:FlowActionPeriodAfter 0.249071
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.374929
## StationCodeI80:Flow:FlowActionPeriodDuring 0.000166 ***
## StationCodeLIS:Flow:FlowActionPeriodDuring 0.231434
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.349200
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.629007
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.868712
## StationCodeRD22:Flow:FlowActionPeriodAfter 0.985860
## StationCodeI80:Flow:FlowActionPeriodAfter 0.002555 **
## StationCodeLIS:Flow:FlowActionPeriodAfter 0.969721
## StationCodeSTTD:Flow:FlowActionPeriodAfter 0.000537 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter 0.762798
## StationCodeRVB:Flow:FlowActionPeriodAfter 0.926785
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.651 3.651 9.374 1.24e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.806
## Scale est. = 0.25102 n = 3478
appraise(m_gamm_q_fa_ar2_gam)
k.check(m_gamm_q_fa_ar2_gam)
## k' edf k-index p-value
## s(DOY) 4 3.650986 0.8149469 0
draw(m_gamm_q_fa_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_fa_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_q_fa_ar2)
Box.test(resid_q_fa_ar2, lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: resid_q_fa_ar2
## X-squared = 30.35, df = 20, p-value = 0.06438
The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_fa_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 0.342 0.849425
## StationCode 6 6.770 3.80e-07
## Flow 1 0.213 0.644739
## FlowActionPeriod 2 0.607 0.544928
## Year_fct:StationCode 24 3.414 3.87e-08
## StationCode:Flow 6 4.543 0.000135
## Flow:FlowActionPeriod 2 0.026 0.974524
## StationCode:FlowActionPeriod 12 5.610 1.29e-09
## StationCode:Flow:FlowActionPeriod 12 5.040 2.22e-08
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.651 3.651 9.374 1.24e-06
Despite its complexity, we’ll use this model for the model selection process.
Finally, we’ll try running a GAM using the model structure from the flow with flow action period model and adding two-way interactions between Year:WaterTemp and Station:WaterTemp. First we’ll run the GAM without accounting for serial autocorrelation.
m_gam_q_fa_wt <- gam(
Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5),
data = df_chla_wt_q,
method = "REML"
)
Lets look at the model summary and diagnostics:
summary(m_gam_q_fa_wt)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 9.547e+00 2.926e-01 32.628
## Year_fct2016 -3.523e-01 3.056e-01 -1.153
## Year_fct2017 -5.764e-03 1.687e-01 -0.034
## Year_fct2018 -3.353e-01 1.608e-01 -2.086
## Year_fct2019 4.571e-01 1.558e-01 2.933
## StationCodeRD22 -4.980e-01 3.355e-01 -1.484
## StationCodeI80 7.357e-01 3.224e-01 2.281
## StationCodeLIS -2.400e+00 2.969e-01 -8.084
## StationCodeSTTD -1.276e+00 3.227e-01 -3.955
## StationCodeLIB -1.247e+00 4.001e-01 -3.117
## StationCodeRVB -2.148e+00 3.514e-01 -6.113
## Flow -7.271e-04 6.433e-04 -1.130
## FlowActionPeriodDuring 3.434e-02 9.046e-02 0.380
## FlowActionPeriodAfter 4.980e-01 1.057e-01 4.710
## WaterTemp -1.568e-02 1.146e-02 -1.368
## Year_fct2016:StationCodeRD22 -4.571e-01 1.005e-01 -4.548
## Year_fct2017:StationCodeRD22 -5.739e-01 8.460e-02 -6.783
## Year_fct2018:StationCodeRD22 -2.862e-01 8.045e-02 -3.558
## Year_fct2019:StationCodeRD22 -6.503e-02 8.367e-02 -0.777
## Year_fct2016:StationCodeI80 3.141e-01 1.040e-01 3.022
## Year_fct2017:StationCodeI80 -1.949e-01 8.623e-02 -2.260
## Year_fct2018:StationCodeI80 -3.712e-01 8.364e-02 -4.438
## Year_fct2019:StationCodeI80 6.428e-01 8.528e-02 7.538
## Year_fct2016:StationCodeLIS 1.408e-01 9.922e-02 1.419
## Year_fct2017:StationCodeLIS 1.859e-01 8.511e-02 2.185
## Year_fct2018:StationCodeLIS 1.916e-02 7.961e-02 0.241
## Year_fct2019:StationCodeLIS 6.092e-01 8.473e-02 7.190
## Year_fct2016:StationCodeSTTD 1.687e-01 1.025e-01 1.647
## Year_fct2017:StationCodeSTTD -5.588e-02 9.667e-02 -0.578
## Year_fct2018:StationCodeSTTD -5.559e-01 8.433e-02 -6.592
## Year_fct2019:StationCodeSTTD -1.018e+00 8.731e-02 -11.662
## Year_fct2016:StationCodeLIB 8.683e-01 1.027e-01 8.458
## Year_fct2017:StationCodeLIB -2.345e-01 9.240e-02 -2.538
## Year_fct2018:StationCodeLIB -1.994e+00 1.094e-01 -18.217
## Year_fct2019:StationCodeLIB -9.703e-01 1.172e-01 -8.279
## Year_fct2016:StationCodeRVB 1.142e-01 1.092e-01 1.046
## Year_fct2017:StationCodeRVB -3.724e-01 9.999e-02 -3.724
## Year_fct2018:StationCodeRVB -2.509e-01 8.781e-02 -2.857
## Year_fct2019:StationCodeRVB -6.336e-01 9.285e-02 -6.824
## StationCodeRD22:Flow 6.978e-04 1.383e-03 0.505
## StationCodeI80:Flow 4.347e-04 1.870e-03 0.232
## StationCodeLIS:Flow 5.075e-03 9.885e-04 5.134
## StationCodeSTTD:Flow 3.462e-03 1.099e-03 3.149
## StationCodeLIB:Flow 1.214e-03 6.495e-04 1.870
## StationCodeRVB:Flow 7.217e-04 6.433e-04 1.122
## Flow:FlowActionPeriodDuring -1.792e-04 6.558e-04 -0.273
## Flow:FlowActionPeriodAfter -6.344e-03 1.178e-03 -5.385
## StationCodeRD22:FlowActionPeriodDuring -1.234e-01 1.156e-01 -1.068
## StationCodeI80:FlowActionPeriodDuring -8.921e-02 1.197e-01 -0.746
## StationCodeLIS:FlowActionPeriodDuring 5.136e-01 1.132e-01 4.539
## StationCodeSTTD:FlowActionPeriodDuring 5.521e-01 1.237e-01 4.463
## StationCodeLIB:FlowActionPeriodDuring -1.717e-01 1.636e-01 -1.050
## StationCodeRVB:FlowActionPeriodDuring -1.081e-01 1.414e-01 -0.764
## StationCodeRD22:FlowActionPeriodAfter -2.807e-01 1.429e-01 -1.965
## StationCodeI80:FlowActionPeriodAfter -5.002e-01 1.462e-01 -3.422
## StationCodeLIS:FlowActionPeriodAfter -4.052e-02 1.220e-01 -0.332
## StationCodeSTTD:FlowActionPeriodAfter -1.510e-01 1.302e-01 -1.160
## StationCodeLIB:FlowActionPeriodAfter -1.600e+00 1.802e-01 -8.878
## StationCodeRVB:FlowActionPeriodAfter -6.906e-01 1.593e-01 -4.336
## Year_fct2016:WaterTemp 1.554e-02 1.201e-02 1.294
## Year_fct2017:WaterTemp 7.768e-03 6.761e-03 1.149
## Year_fct2018:WaterTemp 1.534e-02 6.559e-03 2.339
## Year_fct2019:WaterTemp -2.449e-02 6.218e-03 -3.938
## StationCodeRD22:WaterTemp 4.178e-02 1.249e-02 3.345
## StationCodeI80:WaterTemp -2.552e-03 1.187e-02 -0.215
## StationCodeLIS:WaterTemp 8.755e-02 1.119e-02 7.823
## StationCodeSTTD:WaterTemp 2.874e-02 1.227e-02 2.343
## StationCodeLIB:WaterTemp 2.504e-02 1.611e-02 1.555
## StationCodeRVB:WaterTemp 3.059e-02 1.375e-02 2.224
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.012e-03 1.393e-03 -0.727
## StationCodeI80:Flow:FlowActionPeriodDuring -1.383e-03 1.874e-03 -0.738
## StationCodeLIS:Flow:FlowActionPeriodDuring -5.526e-03 1.007e-03 -5.491
## StationCodeSTTD:Flow:FlowActionPeriodDuring -1.822e-03 1.117e-03 -1.631
## StationCodeLIB:Flow:FlowActionPeriodDuring -1.014e-05 6.641e-04 -0.015
## StationCodeRVB:Flow:FlowActionPeriodDuring 1.844e-04 6.559e-04 0.281
## StationCodeRD22:Flow:FlowActionPeriodAfter -4.200e-03 2.534e-03 -1.657
## StationCodeI80:Flow:FlowActionPeriodAfter -1.946e-03 2.783e-03 -0.699
## StationCodeLIS:Flow:FlowActionPeriodAfter 8.734e-03 1.560e-03 5.598
## StationCodeSTTD:Flow:FlowActionPeriodAfter 1.013e-02 1.662e-03 6.094
## StationCodeLIB:Flow:FlowActionPeriodAfter 5.708e-03 1.184e-03 4.820
## StationCodeRVB:Flow:FlowActionPeriodAfter 6.351e-03 1.178e-03 5.389
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2016 0.249091
## Year_fct2017 0.972744
## Year_fct2018 0.037076 *
## Year_fct2019 0.003375 **
## StationCodeRD22 0.137856
## StationCodeI80 0.022584 *
## StationCodeLIS 8.60e-16 ***
## StationCodeSTTD 7.82e-05 ***
## StationCodeLIB 0.001842 **
## StationCodeRVB 1.09e-09 ***
## Flow 0.258469
## FlowActionPeriodDuring 0.704248
## FlowActionPeriodAfter 2.58e-06 ***
## WaterTemp 0.171506
## Year_fct2016:StationCodeRD22 5.60e-06 ***
## Year_fct2017:StationCodeRD22 1.38e-11 ***
## Year_fct2018:StationCodeRD22 0.000379 ***
## Year_fct2019:StationCodeRD22 0.437081
## Year_fct2016:StationCodeI80 0.002533 **
## Year_fct2017:StationCodeI80 0.023879 *
## Year_fct2018:StationCodeI80 9.35e-06 ***
## Year_fct2019:StationCodeI80 6.10e-14 ***
## Year_fct2016:StationCodeLIS 0.155939
## Year_fct2017:StationCodeLIS 0.028987 *
## Year_fct2018:StationCodeLIS 0.809791
## Year_fct2019:StationCodeLIS 7.93e-13 ***
## Year_fct2016:StationCodeSTTD 0.099692 .
## Year_fct2017:StationCodeSTTD 0.563289
## Year_fct2018:StationCodeSTTD 5.03e-11 ***
## Year_fct2019:StationCodeSTTD < 2e-16 ***
## Year_fct2016:StationCodeLIB < 2e-16 ***
## Year_fct2017:StationCodeLIB 0.011181 *
## Year_fct2018:StationCodeLIB < 2e-16 ***
## Year_fct2019:StationCodeLIB < 2e-16 ***
## Year_fct2016:StationCodeRVB 0.295425
## Year_fct2017:StationCodeRVB 0.000199 ***
## Year_fct2018:StationCodeRVB 0.004298 **
## Year_fct2019:StationCodeRVB 1.05e-11 ***
## StationCodeRD22:Flow 0.613905
## StationCodeI80:Flow 0.816208
## StationCodeLIS:Flow 2.99e-07 ***
## StationCodeSTTD:Flow 0.001652 **
## StationCodeLIB:Flow 0.061637 .
## StationCodeRVB:Flow 0.261987
## Flow:FlowActionPeriodDuring 0.784664
## Flow:FlowActionPeriodAfter 7.73e-08 ***
## StationCodeRD22:FlowActionPeriodDuring 0.285703
## StationCodeI80:FlowActionPeriodDuring 0.455987
## StationCodeLIS:FlowActionPeriodDuring 5.86e-06 ***
## StationCodeSTTD:FlowActionPeriodDuring 8.33e-06 ***
## StationCodeLIB:FlowActionPeriodDuring 0.293788
## StationCodeRVB:FlowActionPeriodDuring 0.444642
## StationCodeRD22:FlowActionPeriodAfter 0.049498 *
## StationCodeI80:FlowActionPeriodAfter 0.000629 ***
## StationCodeLIS:FlowActionPeriodAfter 0.739835
## StationCodeSTTD:FlowActionPeriodAfter 0.246237
## StationCodeLIB:FlowActionPeriodAfter < 2e-16 ***
## StationCodeRVB:FlowActionPeriodAfter 1.49e-05 ***
## Year_fct2016:WaterTemp 0.195712
## Year_fct2017:WaterTemp 0.250655
## Year_fct2018:WaterTemp 0.019415 *
## Year_fct2019:WaterTemp 8.39e-05 ***
## StationCodeRD22:WaterTemp 0.000832 ***
## StationCodeI80:WaterTemp 0.829731
## StationCodeLIS:WaterTemp 6.85e-15 ***
## StationCodeSTTD:WaterTemp 0.019208 *
## StationCodeLIB:WaterTemp 0.120121
## StationCodeRVB:WaterTemp 0.026229 *
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.467471
## StationCodeI80:Flow:FlowActionPeriodDuring 0.460348
## StationCodeLIS:Flow:FlowActionPeriodDuring 4.30e-08 ***
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.103086
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.987816
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.778648
## StationCodeRD22:Flow:FlowActionPeriodAfter 0.097598 .
## StationCodeI80:Flow:FlowActionPeriodAfter 0.484364
## StationCodeLIS:Flow:FlowActionPeriodAfter 2.34e-08 ***
## StationCodeSTTD:Flow:FlowActionPeriodAfter 1.22e-09 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter 1.50e-06 ***
## StationCodeRVB:Flow:FlowActionPeriodAfter 7.55e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.848 3.988 38.09 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.885 Deviance explained = 88.8%
## -REML = 1434.9 Scale est. = 0.11336 n = 3478
appraise(m_gam_q_fa_wt)
k.check(m_gam_q_fa_wt)
## k' edf k-index p-value
## s(DOY) 4 3.84787 0.9605145 0.015
draw(m_gam_q_fa_wt, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gam_q_fa_wt, pages = 1, all.terms = TRUE)
acf(residuals(m_gam_q_fa_wt))
Box.test(residuals(m_gam_q_fa_wt), lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: residuals(m_gam_q_fa_wt)
## X-squared = 6772.7, df = 20, p-value < 2.2e-16
Now, we’ll try to deal with the residual autocorrelation. We’ll run AR(1), AR(2), and AR(3) models and compare them using AIC.
# Define model formula as an object
f_gam_q_fa_wt <- as.formula(
"Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod + Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)"
)
# Fit original model with k = 5 and three successive AR(p) models
m_gamm_q_fa_wt <- gamm(
f_gam_q_fa_wt,
data = df_chla_wt_q,
method = "REML"
)
m_gamm_q_fa_wt_ar1 <- gamm(
f_gam_q_fa_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 1), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_fa_wt_ar2 <- gamm(
f_gam_q_fa_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 2), # grouped by Year_fct and StationCode
method = "REML"
)
m_gamm_q_fa_wt_ar3 <- gamm(
f_gam_q_fa_wt,
data = df_chla_wt_q,
correlation = corARMA(form = ~ 1|Year_fct/StationCode, p = 3), # grouped by Year_fct and StationCode
method = "REML"
)
# Compare models
anova(
m_gamm_q_fa_wt$lme,
m_gamm_q_fa_wt_ar1$lme,
m_gamm_q_fa_wt_ar2$lme,
m_gamm_q_fa_wt_ar3$lme
)
## Model df AIC BIC logLik Test L.Ratio
## m_gamm_q_fa_wt$lme 1 84 3037.848 3552.798 -1434.924
## m_gamm_q_fa_wt_ar1$lme 2 85 -1962.244 -1441.163 1066.122 1 vs 2 5002.092
## m_gamm_q_fa_wt_ar2$lme 3 86 -1985.590 -1458.380 1078.795 2 vs 3 25.347
## m_gamm_q_fa_wt_ar3$lme 4 87 -1983.981 -1450.641 1078.991 3 vs 4 0.391
## p-value
## m_gamm_q_fa_wt$lme
## m_gamm_q_fa_wt_ar1$lme <.0001
## m_gamm_q_fa_wt_ar2$lme <.0001
## m_gamm_q_fa_wt_ar3$lme 0.5318
Again, it looks like the AR(2) model has the best fit according to the AIC values and backed up by the BIC values. Let’s take a closer look at that one.
# Pull out the residuals and the GAM model
resid_q_fa_wt_ar2 <- residuals(m_gamm_q_fa_wt_ar2$lme, type = "normalized")
m_gamm_q_fa_wt_ar2_gam <- m_gamm_q_fa_wt_ar2$gam
summary(m_gamm_q_fa_wt_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric coefficients:
## Estimate Std. Error t value
## (Intercept) 8.796e+00 6.363e-01 13.824
## Year_fct2016 1.238e+00 6.561e-01 1.887
## Year_fct2017 1.791e-01 5.995e-01 0.299
## Year_fct2018 -4.588e-01 6.028e-01 -0.761
## Year_fct2019 -1.297e-01 5.750e-01 -0.226
## StationCodeRD22 -1.322e+00 6.714e-01 -1.969
## StationCodeI80 4.878e-01 6.551e-01 0.745
## StationCodeLIS -1.669e+00 6.586e-01 -2.534
## StationCodeSTTD -3.777e-01 6.810e-01 -0.555
## StationCodeLIB -3.107e-01 7.561e-01 -0.411
## StationCodeRVB -1.787e+00 8.259e-01 -2.164
## Flow -2.257e-04 3.659e-04 -0.617
## FlowActionPeriodDuring 7.161e-02 1.103e-01 0.649
## FlowActionPeriodAfter 1.159e-01 1.356e-01 0.855
## WaterTemp 1.664e-02 2.240e-02 0.743
## Year_fct2016:StationCodeRD22 -6.248e-01 5.824e-01 -1.073
## Year_fct2017:StationCodeRD22 -6.087e-01 5.523e-01 -1.102
## Year_fct2018:StationCodeRD22 -3.091e-01 5.418e-01 -0.571
## Year_fct2019:StationCodeRD22 -1.262e-01 5.437e-01 -0.232
## Year_fct2016:StationCodeI80 1.584e-01 5.929e-01 0.267
## Year_fct2017:StationCodeI80 -2.583e-01 5.628e-01 -0.459
## Year_fct2018:StationCodeI80 -5.353e-01 5.522e-01 -0.969
## Year_fct2019:StationCodeI80 4.972e-01 5.542e-01 0.897
## Year_fct2016:StationCodeLIS 6.262e-01 5.677e-01 1.103
## Year_fct2017:StationCodeLIS 1.909e-01 5.499e-01 0.347
## Year_fct2018:StationCodeLIS -1.089e-01 5.368e-01 -0.203
## Year_fct2019:StationCodeLIS 3.895e-01 5.428e-01 0.718
## Year_fct2016:StationCodeSTTD 3.693e-01 5.864e-01 0.630
## Year_fct2017:StationCodeSTTD -6.079e-01 5.920e-01 -1.027
## Year_fct2018:StationCodeSTTD -9.043e-01 5.621e-01 -1.609
## Year_fct2019:StationCodeSTTD -1.126e+00 5.536e-01 -2.034
## Year_fct2016:StationCodeLIB 9.232e-01 5.719e-01 1.614
## Year_fct2017:StationCodeLIB -3.232e-01 5.559e-01 -0.581
## Year_fct2018:StationCodeLIB -2.101e+00 6.040e-01 -3.478
## Year_fct2019:StationCodeLIB -8.633e-01 6.117e-01 -1.411
## Year_fct2016:StationCodeRVB 2.471e-01 5.819e-01 0.425
## Year_fct2017:StationCodeRVB -4.778e-01 5.461e-01 -0.875
## Year_fct2018:StationCodeRVB -3.624e-01 5.353e-01 -0.677
## Year_fct2019:StationCodeRVB -7.749e-01 5.446e-01 -1.423
## StationCodeRD22:Flow 6.934e-04 1.181e-03 0.587
## StationCodeI80:Flow 4.800e-03 1.591e-03 3.017
## StationCodeLIS:Flow -2.042e-04 9.037e-04 -0.226
## StationCodeSTTD:Flow 1.321e-04 1.033e-03 0.128
## StationCodeLIB:Flow 4.830e-04 3.730e-04 1.295
## StationCodeRVB:Flow 2.179e-04 3.660e-04 0.595
## Flow:FlowActionPeriodDuring 2.747e-05 3.769e-04 0.073
## Flow:FlowActionPeriodAfter 6.981e-05 7.286e-04 0.096
## StationCodeRD22:FlowActionPeriodDuring 7.152e-02 1.515e-01 0.472
## StationCodeI80:FlowActionPeriodDuring 3.482e-01 1.549e-01 2.248
## StationCodeLIS:FlowActionPeriodDuring -7.908e-02 1.334e-01 -0.593
## StationCodeSTTD:FlowActionPeriodDuring 9.413e-03 1.343e-01 0.070
## StationCodeLIB:FlowActionPeriodDuring -4.343e-01 1.791e-01 -2.426
## StationCodeRVB:FlowActionPeriodDuring -1.575e-01 1.846e-01 -0.854
## StationCodeRD22:FlowActionPeriodAfter 2.049e-01 2.013e-01 1.018
## StationCodeI80:FlowActionPeriodAfter 4.531e-01 2.040e-01 2.221
## StationCodeLIS:FlowActionPeriodAfter 1.002e-01 1.688e-01 0.594
## StationCodeSTTD:FlowActionPeriodAfter -4.901e-01 1.705e-01 -2.874
## StationCodeLIB:FlowActionPeriodAfter -4.527e-01 2.014e-01 -2.247
## StationCodeRVB:FlowActionPeriodAfter -1.957e-01 2.062e-01 -0.949
## Year_fct2016:WaterTemp -5.317e-02 2.012e-02 -2.642
## Year_fct2017:WaterTemp 4.584e-03 1.825e-02 0.251
## Year_fct2018:WaterTemp 2.569e-02 1.940e-02 1.325
## Year_fct2019:WaterTemp 8.629e-03 1.730e-02 0.499
## StationCodeRD22:WaterTemp 6.962e-02 2.193e-02 3.175
## StationCodeI80:WaterTemp -8.944e-03 2.043e-02 -0.438
## StationCodeLIS:WaterTemp 4.623e-02 2.165e-02 2.135
## StationCodeSTTD:WaterTemp 5.668e-03 2.295e-02 0.247
## StationCodeLIB:WaterTemp -3.515e-02 2.818e-02 -1.248
## StationCodeRVB:WaterTemp 1.683e-02 3.139e-02 0.536
## StationCodeRD22:Flow:FlowActionPeriodDuring -1.121e-03 1.164e-03 -0.963
## StationCodeI80:Flow:FlowActionPeriodDuring -5.969e-03 1.557e-03 -3.832
## StationCodeLIS:Flow:FlowActionPeriodDuring 4.993e-04 9.440e-04 0.529
## StationCodeSTTD:Flow:FlowActionPeriodDuring 7.451e-04 1.066e-03 0.699
## StationCodeLIB:Flow:FlowActionPeriodDuring -2.590e-04 3.895e-04 -0.665
## StationCodeRVB:Flow:FlowActionPeriodDuring -1.933e-05 3.772e-04 -0.051
## StationCodeRD22:Flow:FlowActionPeriodAfter -9.655e-04 2.553e-03 -0.378
## StationCodeI80:Flow:FlowActionPeriodAfter -7.980e-03 2.766e-03 -2.885
## StationCodeLIS:Flow:FlowActionPeriodAfter -6.062e-04 1.275e-03 -0.475
## StationCodeSTTD:Flow:FlowActionPeriodAfter 4.649e-03 1.391e-03 3.343
## StationCodeLIB:Flow:FlowActionPeriodAfter -1.939e-04 7.341e-04 -0.264
## StationCodeRVB:Flow:FlowActionPeriodAfter -6.629e-05 7.287e-04 -0.091
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## Year_fct2016 0.059184 .
## Year_fct2017 0.765131
## Year_fct2018 0.446654
## Year_fct2019 0.821571
## StationCodeRD22 0.049085 *
## StationCodeI80 0.456585
## StationCodeLIS 0.011330 *
## StationCodeSTTD 0.579142
## StationCodeLIB 0.681142
## StationCodeRVB 0.030529 *
## Flow 0.537306
## FlowActionPeriodDuring 0.516207
## FlowActionPeriodAfter 0.392757
## WaterTemp 0.457427
## Year_fct2016:StationCodeRD22 0.283391
## Year_fct2017:StationCodeRD22 0.270529
## Year_fct2018:StationCodeRD22 0.568362
## Year_fct2019:StationCodeRD22 0.816454
## Year_fct2016:StationCodeI80 0.789413
## Year_fct2017:StationCodeI80 0.646245
## Year_fct2018:StationCodeI80 0.332393
## Year_fct2019:StationCodeI80 0.369731
## Year_fct2016:StationCodeLIS 0.270033
## Year_fct2017:StationCodeLIS 0.728490
## Year_fct2018:StationCodeLIS 0.839244
## Year_fct2019:StationCodeLIS 0.473020
## Year_fct2016:StationCodeSTTD 0.528925
## Year_fct2017:StationCodeSTTD 0.304603
## Year_fct2018:StationCodeSTTD 0.107788
## Year_fct2019:StationCodeSTTD 0.042061 *
## Year_fct2016:StationCodeLIB 0.106552
## Year_fct2017:StationCodeLIB 0.560996
## Year_fct2018:StationCodeLIB 0.000511 ***
## Year_fct2019:StationCodeLIB 0.158224
## Year_fct2016:StationCodeRVB 0.671099
## Year_fct2017:StationCodeRVB 0.381721
## Year_fct2018:StationCodeRVB 0.498377
## Year_fct2019:StationCodeRVB 0.154848
## StationCodeRD22:Flow 0.557060
## StationCodeI80:Flow 0.002576 **
## StationCodeLIS:Flow 0.821271
## StationCodeSTTD:Flow 0.898262
## StationCodeLIB:Flow 0.195464
## StationCodeRVB:Flow 0.551666
## Flow:FlowActionPeriodDuring 0.941913
## Flow:FlowActionPeriodAfter 0.923674
## StationCodeRD22:FlowActionPeriodDuring 0.636993
## StationCodeI80:FlowActionPeriodDuring 0.024636 *
## StationCodeLIS:FlowActionPeriodDuring 0.553410
## StationCodeSTTD:FlowActionPeriodDuring 0.944114
## StationCodeLIB:FlowActionPeriodDuring 0.015339 *
## StationCodeRVB:FlowActionPeriodDuring 0.393441
## StationCodeRD22:FlowActionPeriodAfter 0.308635
## StationCodeI80:FlowActionPeriodAfter 0.026442 *
## StationCodeLIS:FlowActionPeriodAfter 0.552818
## StationCodeSTTD:FlowActionPeriodAfter 0.004072 **
## StationCodeLIB:FlowActionPeriodAfter 0.024681 *
## StationCodeRVB:FlowActionPeriodAfter 0.342716
## Year_fct2016:WaterTemp 0.008277 **
## Year_fct2017:WaterTemp 0.801701
## Year_fct2018:WaterTemp 0.185344
## Year_fct2019:WaterTemp 0.618037
## StationCodeRD22:WaterTemp 0.001511 **
## StationCodeI80:WaterTemp 0.661513
## StationCodeLIS:WaterTemp 0.032853 *
## StationCodeSTTD:WaterTemp 0.804969
## StationCodeLIB:WaterTemp 0.212274
## StationCodeRVB:WaterTemp 0.591847
## StationCodeRD22:Flow:FlowActionPeriodDuring 0.335423
## StationCodeI80:Flow:FlowActionPeriodDuring 0.000129 ***
## StationCodeLIS:Flow:FlowActionPeriodDuring 0.596880
## StationCodeSTTD:Flow:FlowActionPeriodDuring 0.484804
## StationCodeLIB:Flow:FlowActionPeriodDuring 0.506064
## StationCodeRVB:Flow:FlowActionPeriodDuring 0.959130
## StationCodeRD22:Flow:FlowActionPeriodAfter 0.705263
## StationCodeI80:Flow:FlowActionPeriodAfter 0.003938 **
## StationCodeLIS:Flow:FlowActionPeriodAfter 0.634486
## StationCodeSTTD:Flow:FlowActionPeriodAfter 0.000838 ***
## StationCodeLIB:Flow:FlowActionPeriodAfter 0.791721
## StationCodeRVB:Flow:FlowActionPeriodAfter 0.927530
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.623 3.623 7.084 2.26e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.797
## Scale est. = 0.27259 n = 3478
appraise(m_gamm_q_fa_wt_ar2_gam)
k.check(m_gamm_q_fa_wt_ar2_gam)
## k' edf k-index p-value
## s(DOY) 4 3.622926 0.8021592 0
draw(m_gamm_q_fa_wt_ar2_gam, select = 1, residuals = TRUE, rug = FALSE)
plot(m_gamm_q_fa_wt_ar2_gam, pages = 1, all.terms = TRUE)
acf(resid_q_fa_wt_ar2)
Box.test(resid_q_fa_wt_ar2, lag = 20, type = 'Ljung-Box')
##
## Box-Ljung test
##
## data: resid_q_fa_wt_ar2
## X-squared = 30.263, df = 20, p-value = 0.0657
The AR(2) model has much less residual autocorrelation, and the diagnostics plots look pretty good. We’ll take a look at the ANOVA table to see if we need all the interaction terms included in the model.
# the anova.gam function is similar to a type 3 ANOVA
anova(m_gamm_q_fa_wt_ar2_gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Chla_log ~ Year_fct * StationCode + Flow * StationCode * FlowActionPeriod +
## Year_fct * WaterTemp + StationCode * WaterTemp + s(DOY, k = 5)
##
## Parametric Terms:
## df F p-value
## Year_fct 4 2.098 0.078512
## StationCode 6 4.328 0.000234
## Flow 1 0.381 0.537306
## FlowActionPeriod 2 0.366 0.693855
## WaterTemp 1 0.552 0.457427
## Year_fct:StationCode 24 2.863 4.02e-06
## StationCode:Flow 6 3.598 0.001470
## Flow:FlowActionPeriod 2 0.005 0.994855
## StationCode:FlowActionPeriod 12 5.902 2.95e-10
## Year_fct:WaterTemp 4 4.121 0.002471
## StationCode:WaterTemp 6 5.271 2.04e-05
## StationCode:Flow:FlowActionPeriod 12 4.904 4.35e-08
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(DOY) 3.623 3.623 7.084 2.26e-05
Despite its complexity, we’ll use this model for the model selection process.
As a summary, here are the 6 GAM models we are comparing:
m_gamm_cat_ar2 - Categorical predictors only
m_gamm_cat_wt_ar2b - Categorical predictors
adding water temperature m_gamm_q_ar2b - Flow as continuous predictor
m_gamm_q_wt_ar2 - Flow and water temperature
as continuous predictors m_gamm_q_fa_ar2 - Flow with flow action
period m_gamm_q_fa_wt_ar2 - Flow with flow action
period and water temperature # AIC values
df_gamm_aic <-
AIC(
m_gamm_cat_ar2$lme,
m_gamm_cat_wt_ar2b$lme,
m_gamm_q_ar2b$lme,
m_gamm_q_wt_ar2$lme,
m_gamm_q_fa_ar2$lme,
m_gamm_q_fa_wt_ar2$lme
) %>%
as_tibble(rownames = "gamm_model")
# BIC values
df_gamm_bic <-
BIC(
m_gamm_cat_ar2$lme,
m_gamm_cat_wt_ar2b$lme,
m_gamm_q_ar2b$lme,
m_gamm_q_wt_ar2$lme,
m_gamm_q_fa_ar2$lme,
m_gamm_q_fa_wt_ar2$lme
) %>%
as_tibble(rownames = "gamm_model")
# Combine AIC and BIC and display results
left_join(df_gamm_aic, df_gamm_bic, by = join_by(gamm_model, df))
## # A tibble: 6 × 4
## gamm_model df AIC BIC
## <chr> <dbl> <dbl> <dbl>
## 1 m_gamm_cat_ar2$lme 62 -2189. -1809.
## 2 m_gamm_cat_wt_ar2b$lme 65 -2204. -1806.
## 3 m_gamm_q_ar2b$lme 47 -2198. -1909.
## 4 m_gamm_q_wt_ar2$lme 58 -2174. -1818.
## 5 m_gamm_q_fa_ar2$lme 75 -2006. -1546.
## 6 m_gamm_q_fa_wt_ar2$lme 86 -1986. -1458.
According to AIC, model 2 (categorical predictors adding water temperature) was the best model, while BIC preferred model 3 (flow as continuous predictor). We’ll look at these two models more closely.
m_gamm_cat_wt_ar2b %>% saveRDS(here("manuscript_synthesis/model/interim/chla_cat_wt_gamm_model.rds"))
m_gamm_q_ar2b %>% saveRDS(here("manuscript_synthesis/model/interim/chla_flow_gamm_model.rds"))
df_chla_wt_q %>% saveRDS(here("manuscript_synthesis/model/interim/chla_wt_flow_model_data.rds"))